TABLE OF CONTENTS

PHY/riemisot [ Modules ]

[ Top ] [ Modules ]

NOM

 riemisot(rg,ug,phig,alfg,rd,ud,phid,alfd,xi,r,u,p,phi,alf)

DESCRIPTION

 Routine de resolution du probleme de Riemann isotherme
 avec toutes les routines necessaires 
 Routine for solving the isothermal Riemann problem
 issu de V3D (code recherche F77 Helluy/Golay) / from V3D (F77 research code Helluy/Golay)
 la loi de pression est de la forme / the pressure law is of the form
 p=alf p1 + (1-alf) p2
 p1=p0_phy + c_phy**2 (   alf /   phi *rho-rho0A_phy)
 p2=p0_phy + c_phy**2 ((1-alf)/(1-phi)*rho-rho0W_phy)
 alf et phi sont convectes a la vitesse u / alf and phi are convected at velocity u
 
    ENTREES / INPUT
  rg    : masse volumique a gauche / density left
  ug    : vitesse a gauche / velocity left
  phig  : fraction de masse de fluide a gauche / mass fraction of fluid at left
  alfg  : fraction de volume de fluide a gauche / volume fraction of fluid at left
  rd    : masse volumique a droite / density right
  ud    : vitesse a droite / velocity right
  phid  : fraction de masse de fluide a droite / mass fraction of fluid at right
  alfd  : fraction de volume de fluide  a droite / volume fraction of fluid at right
  xi    : x/t   (x/t=0 pour godunov)

    SORTIES / OUTPUT
   r   : masse volumique / density
   u   : vitesse / velocity
   p   : pression / pressure
   phi : fraction de masse de fluide / mass fraction of fluid
   alf : fraction de volume de fluide / volume fraction of fluid

SOURCE

! Ce logiciel est regi par la licence [CeCILL-B]
! This software is governed by the [CeCILL-B] license
!=========================== DEBUT DES DECLARATIONS ==================== 
!.1-----  Implicit, Use
  use phy_typ
  IMPLICIT NONE
!.2-----  Declaration
  real(kind=kind(0.d0)),intent(in)  ::  rg    !! masse volumique a gauche / left density
  real(kind=kind(0.d0)),intent(in)  ::  ug    !! vitesse a gauche / velocity left
  real(kind=kind(0.d0)),intent(in)  ::  phig  !! fraction de masse de fluide a gauche / fraction of fluid mass on left
  real(kind=kind(0.d0)),intent(in)  ::  alfg  !! fraction de volume de fluide a gauche / fraction of fluid volume at left
  real(kind=kind(0.d0)),intent(in)  ::  rd    !! masse volumique a droite / density right
  real(kind=kind(0.d0)),intent(in)  ::  ud    !! vitesse a droite / velocity right
  real(kind=kind(0.d0)),intent(in)  ::  phid  !! fraction de masse de fluide a droite / mass fraction of fluid at right
  real(kind=kind(0.d0)),intent(in)  ::  alfd  !! fraction de volume de fluide  a droite / volume fraction of fluid at right
  real(kind=kind(0.d0)),intent(in)  ::  xi    !! x/t   (x/t=0 pour godunov) / (x/t=0 for godunov)
  real(kind=kind(0.d0)),intent(out) ::  r     !! masse volumique / density
  real(kind=kind(0.d0)),intent(out) ::  u     !! vitesse / velocity
  real(kind=kind(0.d0)),intent(out) ::  p     !! pression / pressure
  real(kind=kind(0.d0)),intent(out) ::  phi   !! fraction de masse de fluide / mass fraction of fluid
  real(kind=kind(0.d0)),intent(out) ::  alf   !! fraction de volume de fluide / volume fraction of fluid 
  
  real(kind=kind(0.d0))             :: r0g,r0d,s1m,s1p,s2m,s2p
  real(kind=kind(0.d0))             :: ugini,udini,um,ui,pmin,pming,pmind
  real(kind=kind(0.d0))             :: pg,pd,pdini,pgini,pm,cd,cg
  real(kind=kind(0.d0))             :: rdini,rgini,r1,r2
  real(kind=kind(0.d0))             :: dp,drg,drd,ff,df
  integer                           :: itermax,iter
  real(kind=kind(0.d0))             :: eps=1.d-10

!=========================== DEBUT DU CODE EXECUTABLE ==================
!-------  Initialisations
  IF (rd.le.0.d0) write(*,*)'rd',rd
  IF (rg.le.0.d0) write(*,*)'rg',rg

  ! vitesses du son à droite et à gauche
  cg=dsqrt(phig*c_phy**2+(1.d0-phig)*c_phy**2)
  cd=dsqrt(phid*c_phy**2+(1.d0-phid)*c_phy**2)

  r0g=(alfg*c_phy**2*rho0A_phy+(1.d0-alfg)*c_phy**2*rho0W_phy)/cg**2
  r0d=(alfd*c_phy**2*rho0A_phy+(1.d0-alfd)*c_phy**2*rho0W_phy)/cd**2

  pming=p0_phy-cg**2*r0g
  pmind=p0_phy-cd**2*r0d
  IF (pming.le.pmind) THEN
     pmin=pmind
  ELSE
     pmin=pming
  END IF

  drg=rg-r0g
  drd=rd-r0d

  pg=p0_phy+cg**2*(rg-r0g)
  pd=p0_phy+cd**2*(rd-r0d)

  p=0.5d0*(pg+pd)

  dp=1
  ff=1
  iter=1
  itermax=100
  
  eps=1.d-6
  
  rgini=rg
  ugini=ug
  pgini=pg
  rdini=rd
  udini=ud
  pdini=pd

  DO WHILE (dabs(dp).gt.eps.and.iter.lt.itermax)

      IF (p.lt.pmin) p=pmin+eps
      iter=iter+1
      r2=(p-p0_phy)/cd**2+r0d
      r1=(p-p0_phy)/cg**2+r0g 
      if(r1.le.0.d0)write(*,*)r1,p,phig
      if(r2.le.0.d0)write(*,*)r2,p,phid
      ff=ud-ug-h(rd,r2,cd)-h(rg,r1,cg)
      df=-dhdr(rd,r2,cd)/cd**2-dhdr(rg,r1,cg)/cg**2
      dp=-ff/df
      p=p+dp

  END DO

  IF (iter.eq.itermax) THEN
      write(*,*) 'Riemiso pas de convergence',dabs(dp/p0_phy)
      write(*,*) 'rg ug phig alfg pg',rg,ug,phig,alfg,pg
      write(*,*) 'rd ud phid alfd pd',rd,ud,phid,alfd,pd
      STOP
  END IF

  ! intermediate densities and velocity
  r2=(p-p0_phy)/cd**2+r0d
  r1=(p-p0_phy)/cg**2+r0g

  um=ug+h(rg,r1,cg)
  um=ud-h(rd,r2,cd)

  ! vitesse d interface
  ui=um
  pm=p

  ! wave velocities

  ! 1-shock or 1-rarefaction
  IF (r1.gt.rg) THEN
      s1m=ug-cg*dsqrt(r1/rg)
      s1p=um-cg*dsqrt(rg/r1)
  ELSE
      s1m=ug-cg
      s1p=um-cg
  END IF

  ! 2-shock or 2-rarefaction
  IF (r2.gt.rd) THEN
      s2m=ud+cd*dsqrt(r2/rd)
      s2p=um+cd*dsqrt(rd/r2)
  ELSE
      s2m=um+cd
      s2p=ud+cd
  END IF

  ! computation of the solution at xi=x/t

  IF (xi.lt.s1m) THEN
      r=rg
      u=ug
      phi=phig
      alf=alfg
      p=pg
  ELSE IF (xi.lt.s1p) THEN
      u=xi+cg
      r=rg*exp((ug-u)/cg)
      phi=phig
      alf=alfg
      p=p0_phy+cg**2*(r-r0g)
  ELSE IF (xi.lt.um) then
      r=r1
      u=um
      phi=phig
      alf=alfg
      p=p0_phy+cg**2*(r-r0g)
  ELSE IF (xi.lt.s2m) then
      r=r2
      u=um
      phi=phid
      alf=alfd
      p=p0_phy+cd**2*(r-r0d)
  ELSE if (xi.lt.s2p) then
      u=xi-cd
      r=rd*exp((u-ud)/cd)
      phi=phid
      alf=alfd
      p=p0_phy+cd**2*(r-r0d)
  ELSE
      r=rd
      u=ud
      phi=phid
      alf=alfd
      p=p0_phy+cd**2*(r-r0d)
  END IF

CONTAINS

!=======================================================================
FUNCTION h(rL,rLocal,c)
!------------------
  use phy
  IMPLICIT NONE

  real(kind=kind(0.d0)), intent(in) :: rL,rLocal,c
  real(kind=kind(0.d0))             :: h

  IF (rLocal*rL.le.0.d0) THEN
      write(*,*)'pb function h: r=',rLocal,' rL=',rL
      STOP
  END IF
  IF (rLocal.gt.rL) THEN
      h=c*(rL-rLocal)/dsqrt(rL*rLocal)
  ELSE
      h=c*log(rL/rLocal)
  END IF
  
END FUNCTION h
!=======================================================================
FUNCTION dhdr(rL,rLocal,c)
!------------------
  use phy
  IMPLICIT NONE

  real(kind=kind(0.d0)), intent(in) :: rL,rLocal,c
  real(kind=kind(0.d0))             :: dhdr,t2,t7
  
  IF(rLocal*rL.le.0.d0) THEN
      write(*,*)'pb function dh'
      STOP
  END IF
  
  IF (rLocal.gt.rL) THEN
      t2 = dsqrt(rL*rLocal)
      t7 = t2*t2
      dhdr = -c/t2-c*(rL-rLocal)/t7/t2*rL/2.D0
  ELSE
      dhdr=-c/rLocal
  END IF
  
END FUNCTION dhdr

END SUBROUTINE riemisot

!===========================   FIN DE LA ROUTINE    ====================