1 : C FLUSTI SOURCE BECC 11/05/26 21:15:25 6981 2 : SUBROUTINE FLUSTI( 3 : & PC_L, GAM_L, 4 : & PC_R, GAM_R, 5 : & Z, 6 : & D_L, D_B, D_A, D_R, 7 : & RHO_L, P_L, U_L, UT1_L, UT2_L, 8 : & RHO_B, P_B, U_B, 9 : & RHO_A, P_A, U_A, 10 : & RHO_R, P_R, U_R, UT1_R, UT2_R, 11 : & UINT, 12 : & FLURHO, FLURU, FLURT1, FLURT2, FLURET, 13 : & LOGDEB, LOGAN) 14 : * 15 : ************************************************************************* 16 : * 17 : * project : CAST3M, EUROPLEXUS... 18 : * 19 : * name : flusti 20 : * 21 : * description : euler equations for a mixture of stiffened gases 22 : * flux in the non-reactive case. 23 : * 24 : * given z = x/t, we compute the flux 25 : * f(x/t) - ((x/t) * u) 26 : * we compute density, pressure and velocity by 27 : * averaging the intermediate states using the 28 : * weve speeds. 29 : * we compute the other variables using these 30 : * quantities and the eos 31 : * 32 : * language : fortran 77 33 : * 34 : * author : a. beccantini den/dm2s/sfme/ltmf 35 : * 36 : ************************************************************************* 37 : * 37 : * called by : 38 : * 39 : * 40 : ************************************************************************* 41 : * 42 : ***** input 43 : * 44 : * pc_l, gam_l = properties of the gas in the left 45 : * 46 : * pc_r, gam_r = properties of the gas in the right 47 : * 48 : * z = x / t 49 : * 50 : * d_l, d_b, d_a, d_r = 51 : * = wave speeds 52 : * 53 : * rho_, p_, u_ 54 : * = density, pressure, velocity 55 : * 56 : * uint = interface speed 57 : * 58 : * logdeb = debugging ? 59 : * 60 : ***** output 61 : * 62 : * flurho, flurhou, flurt1, flurt2, fluret = 63 : * fluxes for rho, rho u, rho ut1, rhout_2, rho e_t 64 : * 65 : * logan = anomaly 66 : * 67 : ************************************************************************* 68 : * 69 : * 26/11/2009 created 70 : * 25/05/2011 evolution in CAST3M 71 : * 72 : ************************************************************************* 73 : * 74 : * n.b.: all variables are declared 75 : * 76 : C IMPLICIT NONE 77 : IMPLICIT INTEGER(I-N) 78 : REAL*8 79 : & PC_L, GAM_L, 80 : & PC_R, GAM_R, 81 : & Z, 82 : & D_L, D_B, D_A, D_R, 83 : & RHO_L, P_L, U_L, UT1_L, UT2_L, 84 : & RHO_B, P_B, U_B, 85 : & RHO_A, P_A, U_A, 86 : & RHO_R, P_R, U_R, UT1_R, UT2_R, 87 : & UINT, 88 : & FLURHO, FLURU, FLURT1, FLURT2, FLURET 89 : * 90 : REAL*8 FP_L, FP_B, FP_A, FP_R 91 : & , GAM_Z, PC_Z 92 : & , RHO_Z, P_Z, U_Z, E_Z, UT1_Z, UT2_Z 93 : & , ECIN_Z 94 : * 95 : * debugging ? 96 : * 97 : LOGICAL LOGDEB, LOGAN 98 : * 99 : ******** states 100 : * 101 : * l 102 : * *** 103 : * |* r 104 : * | * b ********* 105 : * | ******** *| 106 : * | | | * a * | 107 : * | | | ************* | 108 : * | | | | | | 109 : * | | | | | | 110 : * d_l d_b un_b un_a d_a d_r 111 : * 112 : *c 113 : *c test 114 : *c 115 : * z = 1.13d0 116 : * d_l = -1.0d0 117 : * d_b = 1.10d0 118 : * d_a = 1.11d0 119 : * d_r = 1.12d0 120 : * 121 : * speed of the contact discontinuty u_ab = u_a = u_b 122 : * if there is no vacuum between a and b 123 : * 124 : * 125 : ***** weight functions 126 : * 127 : FP_L = 0.0D0 128 : FP_B = 0.0D0 129 : FP_A = 0.0D0 130 : FP_R = 0.0D0 131 : * 132 : * this can be simplified from a computational point of view 133 : * using max, min, abs... 134 : * but in this way it is easy to read it... 135 : * 136 : IF (Z .LT. D_L) THEN 137 : FP_L = 1.0D0 138 : ELSEIF(Z .LT. D_B) THEN 139 : * rarefaction if d_l < d_b 140 : FP_L = (D_B - Z) / (D_B - D_L) 141 : FP_B = 1.0D0 - FP_L 142 : ELSEIF(Z .LT. U_B)THEN 143 : * contact discontinuty, left 144 : FP_B = 1.0D0 145 : ELSEIF(Z .LT. U_A)THEN 146 : * vacuum (u_a < u_b) 147 : FP_B = (U_A - Z) / (U_A - U_B) 148 : FP_A = 1.0D0 - FP_B 149 : ELSEIF(Z .LT. D_A)THEN 150 : * rs, right 151 : FP_A = 1.0D0 152 : ELSEIF(Z .LT. D_R)THEN 153 : * rgnl rarefaction if d_r < d_a 154 : FP_A = (D_R - Z) / (D_R - D_A) 155 : FP_R = 1.0D0 - FP_A 156 : ELSE 157 : * rgnl, right 158 : FP_R = 1.0D0 159 : ENDIF 160 : * 161 : ***** shock-shock + (e.f) intermediate solution 162 : * or 163 : * er-parameterisation, with linearisation over 164 : * the rarefaction for p, rho, w 165 : * 166 : * nb expression of gam_z and pc_z are exact except 167 : * in the case of vacuum (u_b < z < u_a). 168 : * indeed 169 : * gam_z = gam_l if z < un_b 170 : * gam_z = gam_r if z > un_a 171 : * 172 : GAM_Z = ((FP_L + FP_B) * GAM_L) + 173 : & ((FP_A + FP_R) * GAM_R) 174 : PC_Z = ((FP_L + FP_B) * PC_L) + 175 : & ((FP_A + FP_R) * PC_R) 176 : UT1_Z = ((FP_L + FP_B) * UT1_L) + 177 : & ((FP_A + FP_R) * UT1_R) 178 : UT2_Z = ((FP_L + FP_B) * UT2_L) + 179 : & ((FP_A + FP_R) * UT2_R) 180 : RHO_Z = (FP_L * RHO_L) + (FP_B * RHO_B) + 181 : & (FP_A * RHO_A) + (FP_R * RHO_R) 182 : P_Z = (FP_L * P_L) + (FP_B * P_B) + 183 : & (FP_A * P_A) + (FP_R * P_R) 184 : U_Z = (FP_L * U_L) + (FP_B * U_B) + 185 : & (FP_A * U_A) + (FP_R * U_R) 186 : * 187 : IF (LOGDEB) THEN 188 : WRITE(*,*) 'FP_L, FP_B, FP_A, FP_R ', FP_L, FP_B, FP_A, FP_R 189 : WRITE(*,*) 'GAM_Z ', GAM_Z 190 : WRITE(*,*) 'PC_Z ', PC_Z 191 : WRITE(*,*) 'UT1_Z ', UT1_Z 192 : WRITE(*,*) 'UT2_Z ', UT2_Z 193 : WRITE(*,*) 'RHO_Z ', RHO_Z 194 : WRITE(*,*) 'P_Z ', P_Z 195 : WRITE(*,*) 'U_Z ', U_Z 196 : ENDIF 197 : * 198 : ***** computation of the internal energy 199 : * 199 : CALL ESTIFF(RHO_Z, P_Z, GAM_Z, PC_Z, E_Z, LOGDEB, LOGAN) 200 : IF (LOGAN) THEN 201 : WRITE(*,*) 'SUBROUTINE FLUSTI.F' 202 : WRITE(*,*) 'ANOMALY DETECTED' 203 : GOTO 9999 204 : ENDIF 205 : ECIN_Z = 0.5D0 * ((U_Z * U_Z) + (UT1_Z * UT1_Z) + 206 : & (UT2_Z * UT2_Z)) 207 : * 208 : ***** interfacial flux 209 : * according to nkonga, comput methods appl. mech engnr 190, 2000 210 : * \dep{u}{x} + \dep{f(u)}{x} = 0 211 : * z = speed of the surface 212 : * f_z = f - z u 213 : * 214 : UINT = U_Z - Z 215 : FLURHO = RHO_Z * (U_Z - Z) 216 : * write(*,*) 'flurho', flurho 217 : FLURU = (FLURHO * U_Z) + P_Z 218 : FLURET = (FLURHO * (E_Z + ECIN_Z)) 219 : & + (P_Z * U_Z) 220 : FLURT1 = FLURHO * UT1_Z 221 : FLURT2 = FLURHO * UT2_Z 222 : * 223 : * write(*,*) (et_z + ef_z + 224 : * & ecin_z), flurho 225 : * write(*,*) p_z, u_z 226 : 9999 RETURN 227 : END 228 :
© Cast3M 2003 - All rights reserved.
Disclaimer