subroutine SourceTerms_NOTfriction_channel(me,i,QorT)
type(channel_t), intent(inout) :: me
integer, intent(in) :: i
real(dp), intent(in) :: QorT
real(dp) :: eps, ccv, ccp, mmu, llambda, dPdT_Ro, dEdp_Ro, Re
real(dp) :: Vit, Press, e, csound
real(dp) :: dTdp_Ro,dTdRo_p,Der_dTdP_Ro_dRo_T,Der_dTdP_Ro_dT_Ro
real(dp) :: ccv_plusRo,ccp_plusRo,mmu_plusRo,llambda_plusRo,dPdT_Ro_plusRo
real(dp) :: ccv_plusT,ccp_plusT,mmu_plusT,llambda_plusT,dPdT_Ro_plusT
real(dp) :: DerccvdRo_T, DerccvdT_Ro, DerccpdRo_T, DerccpdT_Ro, DermmudRo_T, DermmudT_Ro
real(dp) :: DerllambdadRo_T, DerllambdadT_Ro, Pra, Nu, dNudRe_Pr, dNudPr_Re, dNudPr_Tempk
real(dp) :: Rok,Tempk,Ht,Phik,WetPerim
real(dp), dimension(Nb_VarC) :: dRedUi, Ui, DerpdUi, Der_dEdp_Ro_dUi
real(dp), dimension(Nb_VarC) :: DerccvdUi, DerccpdUi, DermmudUi, DerllambdadUi, Der_dTdp_Ro_dUi
real(dp), dimension(Nb_VarC) :: dPrdUi, DerNudUi
real(dp), dimension(Nb_VarC) :: DerRok, DerHt
real(dp), dimension(Nb_VarC) :: DerPhik, DerTempk
real(dp) :: kk, KKK, dedT, dPdT, HH, dPde
eps=1.0e-2_dp
Ui(:)=me%StVar%He_cs(:,i)
Rok=Ui(Con_Mas)
Vit=Ui(Con_Qdm)/Rok
e=Ui(Con_Ene)/Ui(Con_Mas)-0.5_dp*Vit**2
if(R_Correction) then
call state_roE_withR(Rok, e, Ui(Con_R), Press, Tempk, csound)
else
call state_roE(Rok, e, Press, Tempk, csound)
endif
call he_prop(Rok,Tempk,ccv,ccp,mmu,llambda,dPdT_Ro)
if (sim_error > 0) return
dEdp_Ro=ccv/dPdT_Ro
Re=Rok*abs(Vit)*me%HeProp%Diam/mmu
DerRok(:)=0.0_dp
DerRok(Con_Mas)=1.0_dp
call jacobian_roT(Rok, Tempk, dedT, dPdT, dTdp_Ro, dTdRo_p, &
d2TdP_dT=Der_dTdP_Ro_dT_Ro, d2TdP_dRo=Der_dTdP_Ro_dRo_T)
call he_prop(Rok+eps,Tempk,ccv_plusRo,ccp_plusRo,mmu_plusRo,llambda_plusRo,dPdT_Ro_plusRo)
if (sim_error > 0) return
call he_prop(Rok,Tempk+eps,ccv_plusT,ccp_plusT,mmu_plusT,llambda_plusT,dPdT_Ro_plusT)
if (sim_error > 0) return
DerccvdRo_T=(ccv_plusRo-ccv)/eps
DerccvdT_Ro=(ccv_plusT-ccv)/eps
DerccpdRo_T=(ccp_plusRo-ccp)/eps
DerccpdT_Ro=(ccp_plusT-ccp)/eps
DermmudRo_T=(mmu_plusRo-mmu)/eps
DermmudT_Ro=(mmu_plusT-mmu)/eps
DerllambdadRo_T=(llambda_plusRo-llambda)/eps
DerllambdadT_Ro=(llambda_plusT-llambda)/eps
if(R_Correction) then
DerpdUi(1)=-0.5_dp*(Ui(2)**2)/(Ui(1)**2)
DerpdUi(2)=Ui(2)/Ui(1)
DerpdUi(3)=-1.0_dp
DerpdUi(4)=1.0_dp
else
HH=e+0.5_dp*Vit**2+Press/Rok
call jacobian_roT(Rok, Tempk, dedT, dPdT)
dPde=dPdT/dedT
kk=dPde/Rok
KKK=csound**2+kk*(Vit**2-HH)
DerpdUi(1)=KKK
DerpdUi(2)=-kk*Vit
DerpdUi(3)=kk
DerpdUi(4)=0.0_dp
endif
DerTempk(:)=dTdRo_p*DerRok(:)+dTdp_Ro*DerPdUi(:)
DerccvdUi(:)=DerccvdRo_T*DerRok(:)+DerccvdT_Ro*DerTempk(:)
DerccpdUi(:)=DerccpdRo_T*DerRok(:)+DerccpdT_Ro*DerTempk(:)
DermmudUi(:)=DermmudRo_T*DerRok(:)+DermmudT_Ro*DerTempk(:)
DerllambdadUi(:)=DerllambdadRo_T*DerRok(:)+DerllambdadT_Ro*DerTempk(:)
Der_dTdp_Ro_dUi(:)=Der_dTdP_Ro_dRo_T*DerRok(:)+Der_dTdP_Ro_dT_Ro*DerTempk(:)
Der_dEdp_Ro_dUi(:)=DerccvdUi(:)*dTdp_Ro+ccv*Der_dTdp_Ro_dUi(:)
if(Re>1.0e-8_dp) then
dRedUi(:)=-me%HeProp%Diam*abs(Ui(2))*DermmudUi(:)/(mmu**2)
dRedUi(2)=dRedUi(2)+me%HeProp%Diam*Ui(2)/(mmu*abs(Ui(2)))
else
dRedUi(:)=0.0_dp
endif
Pra=mmu*ccp/llambda
dPrdUi(:)=(llambda*(DermmudUi(:)*ccp+mmu*DerccpdUi(:))-mmu*ccp*DerllambdadUi(:))/(llambda**2)
if(me%HeProp%ExtHeating=='temp') then ! used in general for the heat exchanger treatment
if(associated(me%nuss)) then
Nu=me%nuss%nusselt(Re,Pra,QorT,Tempk)
dNudRe_Pr=me%nuss%nusselt_der_Re(Re,Pra,QorT,Tempk)
dNudPr_Re=me%nuss%nusselt_der_Pra(Re,Pra,QorT,Tempk)
dNudPr_Tempk=me%nuss%nusselt_der_T(Re,Pra,QorT,Tempk)
DerNudUi(:)=dNudRe_Pr*dRedUi(:)+dNudPr_Re*dPrdUi(:)+dNudPr_Tempk*DerTempk(:)
Ht=Nu*llambda/me%HeProp%Diam
DerHt(:)=(1.0_dp/me%HeProp%Diam)*(DerNudUi(:)*llambda+Nu*DerllambdadUi(:))
else
Ht=me%HeProp%Ht
DerHt(:)=0.0_dp
endif
WetPerim=4.0_dp*me%HeProp%Area/me%HeProp%Diam
me%srcPP%Se=WetPerim*Ht*(QorT-Tempk)/me%HeProp%Area
me%srcPP%DerSe(:)=(WetPerim/me%HeProp%Area)*(DerHt(:)*(QorT-Tempk)-Ht*DerTempk(:))
else if(me%HeProp%ExtHeating=='flux') then
me%srcPP%Se=QorT/me%HeProp%Area
me%srcPP%DerSe(:)=0.0_dp
else
me%srcPP%Se=0.0_dp
me%srcPP%DerSe(:)=0.0_dp
endif
Phik=1.0_dp/(Rok*dEdp_Ro)
DerPhik(:)=-(DerRok(:)*dEdp_Ro+Rok*Der_dEdp_Ro_dUi(:))/((Rok*dEdp_Ro)**2)
if(R_Correction) then
me%srcPP%Sr(i)=me%srcPP%Sr(i)+me%srcPP%Se*(1.0_dp+Phik)
me%srcPP%DerSr(:,i)=me%srcPP%DerSr(:,i)+me%srcPP%DerSe(:)*(1.0_dp+Phik)+me%srcPP%Se*DerPhik(:)
endif
if(me%HeProp%ExtHeating=='link') then
! Communication to FS port associated to channel node i
me%HeProp%thermP(i)%p%Pra=Pra
me%HeProp%thermP(i)%p%Re=Re
me%HeProp%thermP(i)%p%dPrdUi(:)=dPrdUi(:)
me%HeProp%thermP(i)%p%dRedUi(:)=dRedUi(:)
me%HeProp%thermP(i)%p%Phi=Phik
me%HeProp%thermP(i)%p%DerPhi(:)=DerPhik(:)
me%HeProp%thermP(i)%p%llambda=llambda
me%HeProp%thermP(i)%p%DerllambdadUi(:)=DerllambdadUi(:)
me%HeProp%thermP(i)%p%TempP=Tempk
me%HeProp%thermP(i)%p%DerTempPdUi(:)=DerTempk(:)
endif
if(me%HeProp%FFsrcLink) then
! Communication to FF_src port associated to channel node i
me%HeProp%thermF(i)%p%Pra=Pra
me%HeProp%thermF(i)%p%Re=Re
me%HeProp%thermF(i)%p%dPrdUi(:)=dPrdUi(:)
me%HeProp%thermF(i)%p%dRedUi(:)=dRedUi(:)
me%HeProp%thermF(i)%p%Phi=Phik
me%HeProp%thermF(i)%p%DerPhi(:)=DerPhik(:)
me%HeProp%thermF(i)%p%llambda=llambda
me%HeProp%thermF(i)%p%DerllambdadUi(:)=DerllambdadUi(:)
me%HeProp%thermF(i)%p%TempP=Tempk
me%HeProp%thermF(i)%p%DerTempPdUi(:)=DerTempk(:)
endif
end subroutine SourceTerms_NOTfriction_channel