phys_prop.f90

Go to the documentation of this file.

module phys_prop
    use constants, only: dp
    implicit none
    private
    integer :: rb_kx=4,rb_iknot=0,rb_inbvx,rb_nx
    real(dp), dimension(:), allocatable :: rb_tx,rb_coef
    integer :: por_kx=2,por_iknot=0,por_inbvx,por_nx
    real(dp), dimension(:), allocatable :: por_tx,por_coef
    integer :: visc_kx=2,visc_iknot=0,visc_inbvx,visc_nx
    real(dp), dimension(:), allocatable :: visc_tx,visc_coef
    public volume_balance,dispress,min_film_thickness,rb_spline_ini,&
        rb,rb_der,porosity_spline_ini,porosity,visc_spline_ini,visc
contains
!********************************BEGINNING*************************************
pure subroutine  volume_balance(y,vt,fs)
    use globals
    real(dp), intent(out) :: fs
    real(dp), intent(out) :: vt
    real(dp), dimension(:), intent(in) :: y
    integer :: i,neq
    real(dp) :: vf,vs
    neq=size(y)
    vf=0
    vs=0
    do i=1,neq
        if (y(i)<strutfilmparameter*y(1)) then
            vf=vf+2*pi*dr*(0.5_dp+i-1)*y(i)*dr
        else
            vs=vs+2*pi*dr*(0.5_dp+i-1)*y(i)*dr
        endif
    enddo
    vs=vs+pi*y(neq)/3*(rd**2+rc*rd+rc**2)-pi*y(neq)*rd**2
    vt=vf+vs
    fs=vs/vt
end subroutine volume_balance
!***********************************END****************************************


!********************************BEGINNING*************************************
pure subroutine dispress(h,dispr,dph)
    use globals
    real(dp), intent(in) :: h
    real(dp), intent(out) :: dispr
    real(dp), intent(out) :: dph
    dispr=bdp*((hdp/h)**(ndp-1)-(hdp/h)**(mdp-1))*(hdp-cdp)
    dph=(bdp*(hdp/h)**mdp*(hdp*mdp+cdp*(h-h*mdp))-&
        bdp*(hdp/h)**ndp*(hdp*ndp+cdp*(h-h*ndp)))/(h*hdp)
end subroutine dispress
!***********************************END****************************************


!********************************BEGINNING*************************************
subroutine min_film_thickness(y,hmin,hloc,havg)
    use globals
    real(dp), dimension(:), intent(in) :: y
    real(dp), intent(out) :: hmin
    real(dp), intent(out) :: hloc
    real(dp), intent(out) :: havg
    integer :: i,n
    hmin=minval(y,dim=1)
    hloc=minloc(y,dim=1)*dr
    havg=0
    n=size(y,1)
    do i=1,n
        if (y(i)<strutfilmparameter*y(1)) then
            havg=havg+2*pi*dr*(0.5_dp+i-1)*y(i)*dr
        else
            exit
        endif
    enddo
    havg=havg/(pi*(dr*(0.5_dp+i-2))**2)
end subroutine min_film_thickness
!***********************************END****************************************


!********************************BEGINNING*************************************
real(dp) function rb_der(t)
    use bspline_module
    real(dp), intent(in) :: t
    integer :: idx,iflag
    idx=1
    call db1val(t,idx,rb_tx,rb_nx,rb_kx,rb_coef,rb_der,iflag,rb_inbvx)
    if (iflag /= 0) then
        print*, 'evaluation of bubble radius derivative from spline failed',&
            iflag
        stop
    endif
endfunction rb_der
!***********************************END****************************************


!********************************BEGINNING*************************************
real(dp) function rb(t)
    use bspline_module
    real(dp), intent(in) :: t
    integer :: idx,iflag
    idx=0
    call db1val(t,idx,rb_tx,rb_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(bblgr_res)
    use bspline_module
    real(dp), dimension(:,:), intent(in) :: bblgr_res
    integer :: iflag
    rb_nx=size(bblgr_res(:,1))
    if (allocated(rb_tx)) deallocate(rb_tx)
    if (allocated(rb_coef)) deallocate(rb_coef)
    allocate(rb_tx(rb_nx+rb_kx),rb_coef(rb_nx))
    rb_tx=0
    call db1ink(bblgr_res(:,1),rb_nx,bblgr_res(:,2),&
        rb_kx,rb_iknot,rb_tx,rb_coef,iflag)
    rb_inbvx=1
    if (iflag /= 0) then
        print*, 'initialization of bubble radius spline failed'
        stop
    endif
end subroutine rb_spline_ini
!***********************************END****************************************


!********************************BEGINNING*************************************
real(dp) function porosity(t)
    use bspline_module
    real(dp), intent(in) :: t
    integer :: idx,iflag
    idx=0
    call db1val(t,idx,por_tx,por_nx,por_kx,por_coef,porosity,iflag,por_inbvx)
    if (iflag /= 0) then
        print*, 'evaluation of porosity from spline failed',iflag
        stop
    endif
endfunction porosity
!***********************************END****************************************


!********************************BEGINNING*************************************
subroutine porosity_spline_ini(bblgr_res)
    use bspline_module
    real(dp), dimension(:,:), intent(in) :: bblgr_res
    integer :: iflag
    por_nx=size(bblgr_res(:,1))
    if (allocated(por_tx)) deallocate(por_tx)
    if (allocated(por_coef)) deallocate(por_coef)
    allocate(por_tx(por_nx+por_kx),por_coef(por_nx))
    por_tx=0
    call db1ink(bblgr_res(:,1),por_nx,bblgr_res(:,3),&
        por_kx,por_iknot,por_tx,por_coef,iflag)
    por_inbvx=1
    if (iflag /= 0) then
        print*, 'initialization of porosity spline failed'
        stop
    endif
end subroutine porosity_spline_ini
!***********************************END****************************************


!********************************BEGINNING*************************************
real(dp) function visc(t)
    use bspline_module
    real(dp), intent(in) :: t
    integer :: idx,iflag
    idx=0
    call db1val(t,idx,visc_tx,visc_nx,visc_kx,visc_coef,visc,iflag,visc_inbvx)
    if (iflag /= 0) then
        print*, 'evaluation of viscosity from spline failed',iflag
        stop
    endif
endfunction visc
!***********************************END****************************************


!********************************BEGINNING*************************************
subroutine visc_spline_ini(bblgr_res)
    use bspline_module
    real(dp), dimension(:,:), intent(in) :: bblgr_res
    integer :: iflag
    visc_nx=size(bblgr_res(:,1))
    if (allocated(visc_tx)) deallocate(visc_tx)
    if (allocated(visc_coef)) deallocate(visc_coef)
    allocate(visc_tx(visc_nx+visc_kx),visc_coef(visc_nx))
    visc_tx=0
    call db1ink(bblgr_res(:,1),visc_nx,bblgr_res(:,4),&
        visc_kx,visc_iknot,visc_tx,visc_coef,iflag)
    visc_inbvx=1
    if (iflag /= 0) then
        print*, 'initialization of viscosity spline failed'
        stop
    endif
end subroutine visc_spline_ini
!***********************************END****************************************
end module phys_prop