junction_dynamic_parameters_fcn Subroutine

public subroutine junction_dynamic_parameters_fcn(me, n, x, fVec)

incoming branch stop index - total number of branches outgoing branch stop index - number of outgoing branches reference branch index incoming branch start index

Type Bound

junction_dynamic_parameters_t

Arguments

Type IntentOptional Attributes Name
class(junction_dynamic_parameters_t), intent(in) :: me
integer, intent(in) :: n

number of incoming branches + 2 * outgoing branches
real(kind=dp), intent(inout) :: x(n)
real(kind=dp), intent(inout) :: fVec(n)

Calls

proc~~junction_dynamic_parameters_fcn~~CallsGraph proc~junction_dynamic_parameters_fcn junction_dynamic_parameters_t%junction_dynamic_parameters_fcn proc~droeint_drop droeint_droP proc~junction_dynamic_parameters_fcn->proc~droeint_drop proc~eos_terms eos_terms proc~droeint_drop->proc~eos_terms proc~t_rop T_roP proc~droeint_drop->proc~t_rop proc~eos_e_terms eos_e_terms proc~eos_terms->proc~eos_e_terms proc~fill_f_terms fill_f_terms proc~eos_terms->proc~fill_f_terms proc~brent brent proc~t_rop->proc~brent proc~t_rop->proc~fill_f_terms proc~fill_g_dreg_terms fill_g_Dreg_terms proc~t_rop->proc~fill_g_dreg_terms f f proc~brent->f proc~set_error set_error proc~brent->proc~set_error proc~zero zero proc~brent->proc~zero proc~fill_e_tpart fill_e_Tpart proc~eos_e_terms->proc~fill_e_tpart proc~fill_ff_terms fill_ff_terms proc~eos_e_terms->proc~fill_ff_terms proc~fill_e_tpart->proc~fill_g_dreg_terms proc~zero->f

Source Code

subroutine junction_dynamic_parameters_fcn(me, n, x, fVec)

    class(junction_dynamic_parameters_t), intent(in) :: me
    integer, intent(in) :: n !! number of incoming branches + 2 * outgoing branches
    real(dp), intent(inout) :: x(n)
    real(dp), intent(inout) :: fVec(n)

    real(dp), dimension(me%NbOut) :: dFdRo, psi, q, Coef
    real(dp) :: EnthalpyMixStar, Denominator
    type(junction_dynamic_star_t) :: star(me%NbTotBr)
    integer :: i, i_out, i_in, i_ref, i_tot, i_kap
    i_tot = me%NbTotBr !! incoming branch stop index - total number of branches
    i_out = me%NbOut !! outgoing branch stop index - number of outgoing branches
    i_ref = i_out + 1  !! reference branch index
    i_in  = i_out + 2  !! incoming branch start index
    i_kap = me%idxKappaRef

    ! Incoming quantities as a function of p(j)
    associate(br => me%br(i_ref:i_tot), st => star(i_ref:i_tot), xx => x(i_ref:i_tot))

        ! Approximate jump relations through the rarefaction wave
        st%vit = br%vit0 - (xx - br%p0 - br%rho0*(br%q0 - br%qBar) + &
                br%dx0*br%Frc/2.0_dp) / (br%rho0 * (br%vit0 - br%SpeedS0))
        st%rho = br%rho0 * (br%vit0 - br%SpeedS0)/(st%vit - br%SpeedS0)

        if(R_Correction) then
            st%rCorr = (br%rCorr0*(br%vit0 - br%SpeedS0) - br%dx0*br%FrcR/2.0_dp) &
                / (st%vit-br%SpeedS0)
            st%ETot = (st%rCorr - xx)/st%rho + 0.5_dp*st%vit**2 + (br%q0 - br%qBar)        
        else
            st%ETot = br%ETot0 + (br%p0*br%vit0 - xx*st%vit) &
                / (br%rho0*(br%vit0 - br%SpeedS0)) + (br%q0-br%qBar)
        endif

        ! Variables involved in the system to solve
        st%q = st%rho * st%vit
        st%Enthalpy = st%ETot + xx/st%rho  ! Total enthalpy

        EnthalpyMixStar = sum(br%AreaJGB * st%q * st%Enthalpy)
        Denominator = sum(br%AreaJGB * st%q)
    end associate

    if(abs(Denominator) < 1.0e-10_dp) Denominator = 1.0e-10_dp
    EnthalpyMixStar = EnthalpyMixStar / Denominator

    ! Outgoing quantities as a function of p(j) and correction term alpha(j)
    ! through the contact discontinuity
    associate(br => me%br(1:i_out), st => star(1:i_out), xx => x(1:i_out))

        ! Approximate jump relations through the shock wave
        st%vit = br%vit0 - (xx - br%p0 - br%rho0*(br%q0 - br%qBar) + &
                br%dx0 * br%Frc / 2.0_dp ) / (br%rho0*(br%vit0 - br%SpeedS0))
        st%rho_int = br%rho0*(br%vit0 - br%SpeedS0) / (st%vit - br%SpeedS0)
        st%ETot_int = br%ETot0 + (br%p0*br%vit0 - xx*st%vit) / (br%rho0 * &
                (br%vit0 - br%SpeedS0)) + (br%q0 - br%qBar)

        st%QDM_int = st%rho_int * st%vit
        st%Energy_int = st%rho_int * st%ETot_int

        ! Correction through the contact discontinuity
        do i = 1, size(dFdRo)
            dFdRo(i) = droeint_droP(st(i)%rho_int, xx(i))
        enddo

        ! Additional term related to the real gas EOS  
        st%rho = st%rho_int + x(i_tot+1:n)
        st%QDM = st%QDM_int + x(i_tot+1:n) * st%vit
        st%Energy = st%Energy_int + x(i_tot+1:n) * (dFdRo + 0.5_dp*st%vit**2)

        ! Variables involved in the system to solve
        st%q = st%QDM
        st%Enthalpy = (st%Energy + xx)/st%rho
    
        ! SYSTEM TO SOLVE : fVec(:)=0.0_dp
        fVec(1) = sum(me%br(1:i_tot)%AreaJGB * star(1:i_tot)%q)

        ! The dat branch is supposed to be the first given incoming pipe
        Coef = 0.0_dp
        where (abs(st%q)>=epsCoef .and. abs(star(i_ref)%q)>=epsCoef)
            q = -br%AreaJGB * st%q / (me%br(i_ref)%AreaJGB * star(i_ref)%q)
            psi = me%br(i_ref)%AreaJGB / br%AreaJGB
            Coef = 1.0_dp - cos((3.0_dp/4.0_dp)*(3.14159_dp-br%ThetaJGB))/(q*psi)
        end where

        fVec(2:i_ref) = x(i_ref) - xx - Coef * st%rho * st%vit**2
        if(me%idxKappaRef == i_ref) &
            fVec(2:i_ref) = x(i_ref) - xx - Coef * me%kappa * st%rho * st%vit**2
        if(me%idxKappaRef < i_ref) &
            fVec(i_kap+1) = x(i_ref) - x(i_kap) - &
                Coef(i_kap) * me%kappa * star(i_kap)%rho * star(i_kap)%vit**2
    end associate
    ! Incoming pipes
    ! Coef_j=0.0_dp ! see Bassett (2003)
    ! - Coef_j*star%rho(i_in:i_tot)*vitStar(i_in:i_tot)**2    
    fVec(i_in:i_tot) = x(i_ref) - x(i_in:i_tot)
    fVec(i_tot+1:n) = star(1:i_out)%Enthalpy - EnthalpyMixStar
end subroutine junction_dynamic_parameters_fcn