init.f90

Go to the documentation of this file.

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

MODULE init

  USE parameters, ONLY : n_vars , n_cry , n_gas , n_mom

  USE parameters, ONLY : idx_p1 , idx_p2 , idx_u1 , idx_u2 , idx_t ,            &
       idx_xd_first , idx_xd_last , idx_alfa_first , idx_alfa_last ,            &
       idx_beta_first , idx_beta_last

  USE parameters, ONLY : idx_mix_mass_eqn , idx_vol1_eqn , idx_mix_mom_eqn ,    &
       idx_rel_vel_eqn , idx_mix_engy_eqn , idx_dis_gas_eqn_first ,             &
       idx_dis_gas_eqn_last , idx_ex_gas_eqn_first , idx_ex_gas_eqn_last ,      &
       idx_cry_eqn_first , idx_cry_eqn_last

  
  IMPLICIT none

  REAL*8 :: alfa1_in
  REAL*8 :: p1_in
  REAL*8 :: p2_in
  REAL*8 :: delta_p_in
  REAL*8 :: u1_in
  REAL*8 :: u2_in
  REAL*8 :: t_in
  REAL*8, ALLOCATABLE :: beta_in(:)
  REAL*8, ALLOCATABLE :: xd_md_in(:)
  REAL*8 :: frag_eff_in

  REAL*8 :: p_out

CONTAINS

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

  SUBROUTINE init_steady(u_0,qp)

    ! external procedures
    USE constitutive, ONLY : f_beta_eq , f_xdis_eq , eval_densities
    USE constitutive, ONLY : f_alfa3, f_alfa

    ! external variables
    USE constitutive, ONLY : p_1 , p_2 , t , x_d_md, x_g, alfa_g_2
    USE constitutive, ONLY : x_ex_dis_in , x_d_md_eq, beta , beta_eq 
    USE constitutive, ONLY : rho_2 , rho_md, rho_g, rho_1, rho_c
    USE constitutive, ONLY : bar_p_m, gamma_m, cv_m
    USE constitutive, ONLY : bar_p_c, gamma_c, cv_c

    IMPLICIT none

    REAL*8, INTENT(IN) :: u_0
    REAL*8, INTENT(OUT) :: qp(n_vars)

    REAL*8 :: r_rho_g(1:n_gas)
    REAL*8 :: r_rho_2
    REAL*8 :: r_rho_1
    REAL*8 :: r_rho_c(1:n_cry)

    REAL*8 :: r_rho_md

    REAL*8 :: r_u_1
    REAL*8 :: r_u_2

    REAL*8 :: r_frag_eff

    REAL*8 :: alfa2_in
    REAL*8 :: alfa_g_in(1:n_gas)
        
    COMPLEX*16 :: x_g_old(1:n_gas)

    REAL*8 :: xd_md_tot

    REAL*8 :: xtot_in

    INTEGER :: i,j,idx,iter,max_iter
    REAL*8 :: error_iter

    p2_in = p1_in + delta_p_in

    ! evaluate the initial crystal volume fraction

    p_1 = dcmplx( p1_in , 0.d0 ) 
    p_2 = dcmplx( p2_in , 0.d0 ) 
    t = dcmplx( t_in , 0.d0 )

    ! evaluate the initial dissolved gas mass fraction
        
    r_rho_md = REAL( ( p_1 + bar_p_m ) / ( T * cv_m * ( gamma_m - 1.D0) ) )
    rho_md = dcmplx(r_rho_md,0.0)
    
    r_rho_c(1:n_cry) = REAL( ( p_1 + bar_p_c(1:n_cry) ) / ( T * cv_c(1:n_cry) * &
         ( gamma_c(1:n_cry) - DCMPLX(1.D0,0.D0) ) ) )    
    
    r_rho_1 = r_rho_md
    
    rho_c(1:n_cry) = dcmplx( r_rho_c(1:n_cry) , 0.d0 )
    rho_1 = dcmplx( r_rho_1 , 0.d0 )
    
    DO i=1,n_gas
       alfa_g_2(i) = dcmplx(1.0 / n_gas,0.0);
       x_g(i) = dcmplx(1e-2,0.0)
    END DO
    
    iter = 1
    max_iter =1000
    error_iter = 1.0
    
    DO WHILE( error_iter .GT. 1e-15 .AND. (iter .LT. max_iter) )
       
       x_g_old = x_g
       
       CALL f_xdis_eq
       
       DO i=1,n_gas
          
          xd_md_in(i) = min( REAL(x_d_md_eq(i)) , x_ex_dis_in(i) )
          
       END DO
       
       ! required by eval_densities
       x_d_md(1:n_gas) = dcmplx( xd_md_in(1:n_gas) , 0.d0 )
       
       CALL f_beta_eq

       beta_in(1:n_cry) = REAL(beta_eq(1:n_cry))
       
       r_u_1 = u_0
       r_u_2 = u_0 + 1.d-10
       
       ! evaluate the volume fractions of the two phases
       
       beta = dcmplx( beta_in(1:n_cry) , 0.d0 )
       
       CALL eval_densities
       
       r_rho_1 = REAL( rho_1 )
       r_rho_2 = REAL( rho_2 )
       r_rho_g = REAL( rho_g )
       
       r_rho_md = REAL( rho_md )
       
       xd_md_tot = sum( xd_md_in(1:n_gas) )
       
       xtot_in = sum( x_ex_dis_in(1:n_gas) ) 
       
       IF ( n_gas .EQ. 1 ) THEN
          
          CALL f_alfa( x_ex_dis_in(1:n_gas) , xd_md_in(1:n_gas) ,               &
               beta_in(1:n_cry) , r_rho_md , r_rho_2 , alfa_g_in(1:n_gas) )
          
       ELSE
          
          CALL f_alfa3( p2_in, x_ex_dis_in(1:n_gas) , beta_in(1:n_cry) ,  &
               r_rho_md , r_rho_g(1:n_gas) , alfa_g_in(1:n_gas) )
          
       END IF
       
       
       DO i = 1,n_gas
          
          alfa_g_in(i) = max( alfa_g_in(i) , 1.d-10 )
          
       END DO
       
       alfa2_in = sum( alfa_g_in(1:n_gas) )
       alfa1_in = 1.d0 - alfa2_in
       
       x_g = alfa_g_in * r_rho_g / ( alfa1_in * r_rho_1 + alfa2_in * r_rho_2 )
       
       alfa_g_2 = alfa_g_in / alfa2_in
       
       error_iter = (maxval(abs(x_g - x_g_old)))
       
    END DO
    
    IF (iter .GE. max_iter) THEN
       
       WRITE(*,*) 'No convergence for initial gas dissolve mass fraction'
       stop
       
    END IF
    
    ! initialize the fragmentation efficiency to zero
    
    r_frag_eff = frag_eff_in

    ! -------- define the indexes of variables and equations --------------------
    
    idx_p1 = 1
    idx_p2 = 2
    idx_u1 = 3
    idx_u2 = 4
    idx_t = 5
    idx_xd_first = 5+1
    idx_xd_last = 5+n_gas
    idx_alfa_first = 5+n_gas+1
    idx_alfa_last = 5+2*n_gas
    idx_beta_first = 5+2*n_gas+1
    idx_beta_last = 5+2*n_gas+n_cry*n_mom

    idx_mix_mass_eqn = 1
    idx_vol1_eqn = 2
    idx_mix_mom_eqn = 3
    idx_rel_vel_eqn = 4
    idx_mix_engy_eqn = 5
    idx_dis_gas_eqn_first = 5+1
    idx_dis_gas_eqn_last = 5+n_gas
    idx_ex_gas_eqn_first = 5+n_gas+1
    idx_ex_gas_eqn_last = 5+2*n_gas 
    idx_cry_eqn_first = 5+2*n_gas+1
    idx_cry_eqn_last = 5+2*n_gas+n_cry*n_mom

    ! --------- define the vector of primitive variables ------------------------
    
    ! Frist phase pressure
    qp(idx_p1) = p1_in
    
    ! Second phase pressure
    qp(idx_p2) = p2_in

    ! Frist phase velocity
    qp(idx_u1) = r_u_1

    ! Second phase velocity
    qp(idx_u2) = r_u_2

    ! Temperature
    qp(idx_t) = t_in

    ! Dissolved gas mass fractions
    idx = idx_xd_first
    DO i = 1,n_gas

       qp(idx) = xd_md_in(i)

       idx = idx + 1

    END DO

    ! Exsolved gas volume fraction
    idx = idx_alfa_first
    DO i = 1,n_gas
              
       qp(idx) = alfa_g_in(i)

       idx = idx + 1

    END DO
    
    ! Crystal volume fractions
    idx = idx_beta_first

    IF ( n_mom .LE. 1 ) THEN

       DO i = 1,n_cry

          qp(idx) = beta_in(i)
          
          idx = idx + 1
          
       END DO

    ELSE

       DO i = 1,n_cry

          DO j = 0,n_mom-1
             
             ! qp(idx) = mom_cry_in(i,j)
             
             idx = idx + 1
       
          END DO

       END DO

    END IF
       
  END SUBROUTINE init_steady

END MODULE init

---

Generated on Tue Mar 7 2017 09:55:58 for MAMMA by   1.8.11