solver_flow.f90

Go to the documentation of this file.

!********************************************************************************
!********************************************************************************

MODULE solver_module
  USE particles_module, ONLY: iclass
  USE current_module
  USE mixture_module
  USE envi_module, ONLY: p, t_a, gi, gas_constair, c_vair

  !
  IMPLICIT NONE
  !
  ! ... Right-Hand Side (rhs) and Known Term (f) of the equation system
  !
  REAL*8, ALLOCATABLE, DIMENSION(:) :: rhstemp
  REAL*8, ALLOCATABLE, DIMENSION(:) :: rhs
  REAL*8, ALLOCATABLE, DIMENSION(:) :: ftemp
  REAL*8, ALLOCATABLE, DIMENSION(:) :: f
  REAL*8, ALLOCATABLE, DIMENSION(:) :: f_stepold
  REAL*8, ALLOCATABLE, DIMENSION(:) :: f_oldold

  REAL*8, ALLOCATABLE, DIMENSION(:) :: f_new
  REAL*8, ALLOCATABLE, DIMENSION(:) :: f_RK
  REAL*8, ALLOCATABLE, DIMENSION(:,:) :: rhs_RK
  REAL*8, ALLOCATABLE, DIMENSION(:) :: rhs4
  REAL*8, ALLOCATABLE, DIMENSION(:) :: f_err

  REAL*8 :: eps_rel , eps_abs

  INTEGER :: itotal

  INTEGER :: n_RK

  REAL*8 :: dh_dr
  REAL*8 :: dbeta_dr
  REAL*8 :: dr
  REAL*8 :: dr0
  REAL*8 :: dr_old
  !
  LOGICAL :: vel_equation

  SAVE

CONTAINS

  !******************************************************************************
  !******************************************************************************

  SUBROUTINE allocate_matrix

    IMPLICIT NONE

    itotal = iclass + 5

    n_rk = 6

    ALLOCATE(rhstemp(itotal))
    ALLOCATE(rhs(itotal))
    ALLOCATE(ftemp(itotal)) 
    ALLOCATE(f(itotal)) 
    ALLOCATE(f_stepold(itotal))   

    ALLOCATE(f_new(itotal))   
    ALLOCATE(f_rk(itotal))   
    ALLOCATE(rhs_rk(n_rk,itotal))
    ALLOCATE(rhs4(itotal))
    ALLOCATE(f_err(itotal))

    f = 0.d0
    ftemp = 0.d0
    rhs = 0.d0
    rhstemp = 0.d0

    RETURN
  END SUBROUTINE allocate_matrix

  !*****************************************************************************
  !*****************************************************************************

  SUBROUTINE rate(rhs_)

    ! External variables
    USE envi_module, ONLY: alpha, theta, fric
    USE particles_module, ONLY: rhosol , diam, solidmassflux_fract, v_s, s,     &
         sumsed, fracsolid, c_d

    IMPLICIT NONE

    REAL*8, DIMENSION(:), INTENT(OUT) :: rhs_

    REAL*8 :: var , coeff_radial
    
    REAL*8 :: term1 , term2

    REAL*8 :: rhs_2a , rhs_2b


    LOGICAL :: no_entr , const_volume_flux

    INTEGER :: i

    IF ( oned_model ) THEN

       var = 1.d0

    ELSE

       var = r

    END IF

    solid_mass_flux = (beta * u * h * var) * ( 1.d0 - n )  
    ! provides the mass_flux of a particle of a given size
    solidmassflux_fract(1:iclass) = solid_mass_flux*fracsolid(1:iclass)  

    ! SEDIMENTATION for the particle fraction
    v_s(1:iclass) = dsqrt( (rhosol(1:iclass) * gi * diam(1:iclass)) /           &
         (c_d(1:iclass) * beta) )

    s(1:iclass) = 1.d0 / (h * u) * solidmassflux_fract(1:iclass) * v_s(1:iclass)

    sumsed = sum(s(1:iclass)) ! sum of sedimentation across all particles sizes
    
    !---- Mass conservation of the mixture   (Eq.1 Bursik & Woods 1996)

    entrainment_rate = epsilon * alpha * u * var
    rhs_(1) = entrainment_rate - sumsed
    
    !---- Momentum Conservation (Eq.4 Bursik & Woods 1996)

    IF ( r_new .NE. r_old ) THEN

       dbeta_dr = (beta_new - beta_old) / (r_new - r_old)
       dh_dr = (h_new - h_old) / (r_new - r_old)

    ELSE

       dbeta_dr = 0.d0
       dh_dr = 0.d0

    END IF

    coeff_radial = ( var - 1.d0 ) / ( r - 1.d0 )

    no_entr = .false.
    
    const_volume_flux = .false.

    IF ( no_entr ) THEN

       term1 = 0.d0

    ELSE
       
       term1 = - u * epsilon * alpha / ( h * beta )

    END IF

    IF ( const_volume_flux ) THEN

       term2 = 0.d0

    ELSE

       term2 = - ( epsilon * alpha * u * var - sumsed ) / ( beta * h * var ) *  &
            dcos(theta) + u / beta * dcos(theta) * dbeta_dr

    END IF

    ! RHS term for the modified momentum equation 
    rhs_2a = ( term1 + ri * ( term2 + coeff_radial * u / var * dcos(theta) +    &
         u / h * dsin(theta) - u / ( 2.d0 * ( beta - alpha ) ) * dcos(theta) *  &
         dbeta_dr ) - fric * u / h ) / ( 1.d0 - ri * cos(theta) )

    ! RHS term for the original equation
    rhs_2b =  - epsilon * alpha * u / ( beta * h ) + ri * u / h * ( - dh_dr     &
         * dcos(theta) + dsin(theta) ) - ri * u  / ( 2.d0 * ( beta - alpha ))   &
         * dcos(theta) * dbeta_dr - fric * u / h

    IF ( vel_equation ) THEN

       rhs_(2) = rhs_2a
       
    ELSE

       rhs_(2) = rhs_2b

    END IF

    !---- Thermal Energy Conservation (Eq.5 Bursik & Woods 1996)

    rhs_(3) = epsilon * alpha * u * var * (c_vair * t_a + p/alpha) - sumsed *   &
         c_s * t    

    !---- Mass average specific heat
    rhs_(4) = ((c_vair - c_vmix) / beta * epsilon * alpha * u * var + (c_s -    &
         c_vmix) / beta * (-sumsed))/ (u * h * var)    

    !---- Gas constant of mixture
    rhs_(5) = (gas_constair - gas_constmix) / (beta * n * u * h * var ) *       &
         epsilon * alpha * u * var

    DO i=1, iclass

       rhs_(5+i) = -s(i)

    ENDDO

    RETURN
  END SUBROUTINE rate


  !*****************************************************************************
  !*****************************************************************************

  SUBROUTINE lump(f_)

    ! External variables
    USE envi_module, ONLY: p
    USE mixture_module, ONLY: c_vmix, cpwvapour
    USE particles_module, ONLY: fracsolid, iclass

    IMPLICIT NONE
    REAL*8, DIMENSION(:), INTENT(OUT) :: f_
    INTEGER :: i
    REAL*8 :: var

    CALL compute_mixture
    !

    IF ( oned_model ) THEN

       var = 1.d0

    ELSE

       var = r

    END IF


    f_(1) = beta * u * h * var

    f_(2) = u
    f_(3) = beta * u * h * var * (c_vmix * t + p/beta + 0.5 * u**2)
    f_(4) = c_vmix
    f_(5) = gas_constmix

    DO i=1, iclass  
     
       f_(5+i) = (beta * u * h * var) * (1-n) * fracsolid(i) 
   
    ENDDO

    RETURN
  END SUBROUTINE lump

  !******************************************************************************
  !
  !*****************************************************************************

  SUBROUTINE marching(fold,fnew,rat,delta_r)

    IMPLICIT NONE
    REAL*8, DIMENSION(:), INTENT(IN) :: fold, rat
    REAL*8, INTENT(IN) :: delta_r
    REAL*8, DIMENSION(:), INTENT(OUT) :: fnew
    INTEGER :: i

    DO i=1,itotal

       fnew(i) = fold(i) + rat(i) * delta_r

    END DO

    RETURN

  END SUBROUTINE marching

  !*****************************************************************************
  !*****************************************************************************

  SUBROUTINE unlump(f_)

    ! External varuables
    USE envi_module, ONLY: alpha , p
    USE mixture_module, ONLY: rhosol_ave
    USE particles_module, ONLY: iclass, fracsolid, rhosol

    IMPLICIT NONE

    REAL*8, DIMENSION(:), INTENT(IN) :: f_

    REAL*8 :: c1 , c2 , c3 , c4

    REAL*8 :: var

    IF ( oned_model ) THEN

       var = 1.d0

    ELSE

       var = r

    END IF

    u = f_(2)
    n = 1.d0 - sum(f_(5+1:5+iclass))/f_(1)    
    c_vmix = f_(4)    
    gas_constmix = f_(5)
    fracsolid(1:iclass) = (f_(5+1:5+iclass))/sum(f_(5+1:5+iclass))

    rhosol_ave = 1.d0 / (sum(fracsolid/rhosol))

    c1 = ( 1.d0 - n ) / rhosol_ave
    c2 = n * gas_constmix / p
    c3 = ( f_(3)/f_(1) - 0.5d0 * u**2 ) / c_vmix
    c4 = - p / c_vmix
    
    beta = ( 1.d0 - c2*c4 ) / ( c2*c3 + c1 )
    t = c4 / beta + c3

    h = f_(1) / (beta * u * var)

    ri = ((beta - alpha) * gi * h) / (beta * u**2.d0) 
    epsilon = 0.075d0 / dsqrt(1.d0 + (718.d0 * ri**2.4d0))

    RETURN

  END SUBROUTINE unlump

END MODULE solver_module