Initial commit from Polytechnique Montreal

1 files changed, 100 insertions, 0 deletions

diff --git a/Donjon/src/THMPH.f b/Donjon/src/THMPH.f
new file mode 100644
index 0000000..487262e
--- /dev/null
+++ b/Donjon/src/THMPH.f
@@ -0,0 +1,100 @@
+*DECK THMPH
+      SUBROUTINE THMPH(IFLUID,PP,HH,RHO,TEMP)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Backwards inversion of steam tables to find water density and
+* temperature.
+*
+*Copyright:
+* Copyright (C) 2012 Ecole Polytechnique de Montreal.
+*
+*Author(s): 
+* A. Hebert, P. Gallet
+*
+*Parameters: input
+* IFLUID  type of fluid (0=H2O; 1=D2O).
+* PP      pressure (Pa)
+* HH      enthalpy (J/Kg)
+*
+*Parameters: output
+* RHO     water density (Kg/m^3)
+* TEMP    temperature (K)
+*
+*-----------------------------------------------------------------------
+*
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      INTEGER IFLUID
+      REAL PP,HH,RHO,TEMP
+*----
+*  LOCAL VARIABLES
+*----
+      PARAMETER(DT=-0.01,ZKELV=273.15,S=1.0)
+      REAL A(15),R1,R3,R4,R5,RV,HSAT,HV,XTH
+      DATA A/
+     > 0.2873E+03,-0.5098E+00,-0.3459E+00,0.1910E+00,-0.2840E-01,
+     > 0.8266E+02,0.1141E+01,-0.2724E+01,0.1077E+00,-0.1144E+02,
+     > 0.9500E+01,-0.2715E+01,-0.1290E+02,0.9148E+01,-0.8093E+01 /
+*----
+*  INITIAL APPROXIMATION OF T1
+*----
+      IF(IFLUID.EQ.0) THEN
+         CALL THMSAT(PP,TSAT)
+         CALL THMTX(TSAT,1.0,RV,HV,R3,R4,R5)
+         CALL THMTX(TSAT,0.0,R1,HSAT,R3,R4,R5)
+      ELSE IF(IFLUID.EQ.1) THEN
+         CALL THMHST(PP,TSAT)
+         CALL THMHTX(TSAT,1.0,RV,HV,R3,R4,R5)
+         CALL THMHTX(TSAT,0.0,R1,HSAT,R3,R4,R5)
+      ENDIF
+      IF((ABS(HSAT-HH)/HSAT).LE.1.0E-5) THEN
+         T1=TSAT
+         GO TO 20
+      ELSEIF(HH.LE.HSAT) THEN
+         T=(HH-1270.0E3)/420.0E3
+         P=(PP-140.0E5)/30.0E5
+         H1=A(1)+P*(A(2)+P*(A(3)+P*(A(4)+P*A(5))))
+         H2=A(6)+P*(A(7)+P*(A(8)+P*(A(9))))
+         H3=A(10)+P*(A(11)+P*(A(12)))
+         H4=A(13)+P*A(14)
+         H5=A(15)
+         T1=H1+T*(H2+T*(H3+T*(H4+T*(H5))))+ZKELV
+* INLET TEMPERATURE WAS VERIFIED TO BE GREATER THAN 0 C. T1 INITIAL
+* GUESS LOWER THAN THAT SHOULD BE INTERPRETED AS FLAWED (FAR FROM
+* FITTING REGION). CORRECTING WITH AN ABOVE-0 C GUESS.
+         IF(T1.LT.ZKELV) T1=10.0+ZKELV
+      ELSEIF(HH.LE.HV) THEN
+*        saturated steam
+         TEMP=TSAT
+         XTH=(HH-HSAT)/(HV-HSAT)
+         RHO=1.0/(XTH/RV+(1.0-XTH)/R1)
+         GO TO 30
+      ELSE
+*        superheated steam
+         T1=TSAT
+      ENDIF
+*----
+*  NEWTON ITERATIONS
+*----
+      ITER=0
+   10 ITER=ITER+1
+      IF(ITER.GT.30) CALL XABORT('THMPH: CONVERGENCE FAILURE.')
+      IF(IFLUID.EQ.0) THEN
+         CALL THMPT(PP,T1,R1,H1,R3,R4,R5)
+         CALL THMPT(PP,T1+DT,R1,H1P,R3,R4,R5)
+      ELSE IF(IFLUID.EQ.1) THEN
+         CALL THMHPT(PP,T1,R1,H1,R3,R4,R5)
+         CALL THMHPT(PP,T1+DT,R1,H1P,R3,R4,R5)
+      ENDIF
+      IF(ABS((HH-H1)/HH).LT.1.E-05) GO TO 20
+      T1=T1+(HH-H1)*DT/(H1P-H1)
+      IF((HH.LE.HSAT).AND.(T1.GE.TSAT)) T1=TSAT
+      IF((HH.GE.HV).AND.(T1.LE.TSAT)) T1=TSAT
+      GO TO 10
+   20 RHO=R1
+      TEMP=T1
+   30 RETURN
+      END