1 : C FLUFU1 SOURCE BECC 09/12/07 21:15:21 6579 2 : subroutine flufu1(nordpo, 3 : & Tmaxcv, Rgas_r, acv_r, Rgas_b, acv_b, 4 : & Rgas_l, acv_l, acv_z, 5 : & z, 6 : & d_l, d_d, d_c, d_b, d_a, d_r, 7 : & rho_l, T_l, u_l, hf_l, 8 : & rho_d, T_d, u_d, 9 : & rho_c, T_c, u_c, 10 : & rho_b, T_b, u_b, hf_b, 11 : & rho_a, T_a, u_a, 12 : & rho_r, T_r, u_r, hf_r, 13 : & flurho, fluru, fluret) 14 : C 15 : c************************************************************************ 16 : c 17 : c project : 18 : c 19 : c name : flufu1 20 : c 21 : c description : Reactive Euler equations for a mixture of reactive 22 : C thermally perfect gases. 23 : c 24 : c Given z = x/t, we compute the flux 25 : C F(x/t) - ((x/t) * U) 26 : C We compute density, temperature and velocity by 27 : C averaging the intermediate states using the 28 : C weve speeds. 29 : C We compute the other variables using these 30 : C quantities and the EOS 31 : c 32 : c language : fortran 77 33 : c 34 : c author : a. beccantini den/dm2s/sfme/ltmf 35 : c 36 : c************************************************************************ 37 : c 37 : C called by : fbecre 38 : C 39 : c 40 : c************************************************************************ 41 : c 42 : c**** input 43 : c 44 : c nordpo = polynomial order of the cv 45 : c 46 : C Tmaxcv = Tmax for cv expansion 47 : C 48 : C z = x / t 49 : C 50 : C Tmaxcv, Rgas_r, acv_r, Rgas_b, acv_b, 51 : C Rgas_l, acv_l = quantities to compute cv and cp 52 : C 53 : C acv_z = work vector 54 : C 55 : C d_l, d_d, un_d, un_c, d_c, d_b, d_a, d_r = 56 : C = wave speeds 57 : C 58 : C rho_, u_, T_, hf_ 59 : C = density, velocity, temperature, 60 : C formation enthalpy 61 : C**** output 62 : C 63 : c flurho, flurhou, fluret = fluxes for rho, rho u, rho e_t 64 : C 65 : c 66 : c************************************************************************ 67 : c 68 : C 07/12/2009 created 69 : C 70 : c************************************************************************ 71 : c 72 : c n.b.: all variables are declared 73 : c 74 : c implicit none 75 : integer nordpo, iord 76 : real*8 Tmaxcv, Rgas_r, acv_r(0:nordpo) 77 : & , Rgas_b, acv_b(0:nordpo) 78 : & , Rgas_l, acv_l(0:nordpo) 79 : & , z, acv_z(0:nordpo) 80 : & , d_l, d_d, d_c, d_b, d_a, d_r 81 : & , rho_l, u_l, T_l, hf_l 82 : & , rho_d, T_d, u_d 83 : & , rho_c, T_c, u_c 84 : & , rho_b, T_b, u_b, hf_b 85 : & , rho_a, T_a, u_a 86 : & , rho_r, T_r, u_r, hf_r 87 : & , flurho, fluru, fluret 88 : real*8 u_cd, fp_l, fp_d, fp_c, fp_b, fp_a, fp_r 89 : & , Rgas_z, rho_z, T_z, u_z, p_z, e_z, cv_z, hf_z 90 : & , ecin_z 91 : logical logdeb 92 : parameter(logdeb = .false.) 93 : CC 94 : CC Test 95 : CC 96 : C z = 1.13D0 97 : C d_l = -1.0D0 98 : C d_d = 0.10D0 99 : C u_c = 0.2D0 100 : C u_d = 0.2D0 101 : C d_c = 1.09D0 102 : C d_b = 1.10D0 103 : C d_a = 1.11D0 104 : C d_r = 1.12D0 105 : C 106 : C Speed of the contact discontinuty u_cd = u_c = u_d 107 : C 108 : u_cd = 0.5D0 * (u_c + u_d) 109 : C 110 : C**** Weight functions 111 : C 112 : fp_l = 0.0d0 113 : fp_d = 0.0D0 114 : fp_c = 0.0d0 115 : fp_b = 0.0d0 116 : fp_a = 0.0d0 117 : fp_r = 0.0d0 118 : C 119 : C This can be simplified from a computational point of view 120 : C 121 : if (z .lt. d_l) then 122 : fp_l = 1.0d0 123 : elseif(z .lt. d_d) then 124 : C Rarefaction if d_l < d_d 125 : fp_l = (d_d - z) / (d_d - d_l) 126 : fp_d = 1.0d0 - fp_l 127 : elseif(z .lt. u_cd)then 128 : C Contact discontinuty, left 129 : fp_d = 1.0d0 130 : elseif(z .lt. d_c)then 131 : C Contact discontinuty, right 132 : fp_c = 1.0d0 133 : elseif(z .lt. d_b)then 134 : C Rarefaction if d_b < d_c 135 : fp_c = (d_b - z) / (d_b - d_c) 136 : fp_b = 1.0d0 - fp_c 137 : elseif(z .lt. d_a)then 138 : C RS, right 139 : fp_a = 1.0D0 140 : elseif(z .lt. d_r)then 141 : C Rgnl rarefaction if d_r < d_a 142 : fp_a = (d_r - z) / (d_r - d_a) 143 : fp_r = 1.0d0 - fp_a 144 : else 145 : C Rgnl, right 146 : fp_r = 1.0d0 147 : endif 148 : CC 149 : C 150 : C**** Shock-shock + (e.f) intermediate solution 151 : C 152 : Rgas_z = ((fp_l + fp_d) * Rgas_l) + 153 : & ((fp_c + fp_b) * Rgas_b) + 154 : & ((fp_a + fp_r) * Rgas_r) 155 : rho_z = (fp_l * rho_l) + (fp_d * rho_d) + 156 : & (fp_c * rho_c) + (fp_b * rho_b) + 157 : & (fp_a * rho_a) + (fp_r * rho_r) 158 : T_z = (fp_l * T_l) + (fp_d * T_d) + 159 : & (fp_c * T_c) + (fp_b * T_b) + 160 : & (fp_a * T_a) + (fp_r * T_r) 161 : p_z = Rgas_z * rho_z * T_z 162 : u_z = (fp_l * u_l) + (fp_d * u_d) + 163 : & (fp_c * u_c) + (fp_b * u_b) + 164 : & (fp_a * u_a) + (fp_r * u_r) 165 : if (logdeb) then 166 : write(*,*) 'd_l, d_d, d_c, d_b, d_a, d_r' 167 : write(*,'(6E12.4)') d_l, d_d, d_c, d_b, d_a, d_r 168 : write(*,*) 'fp_l, fp_d, fp_c, fp_b, fp_a, fp_r' 169 : write(*,'(6E12.4)') fp_l, fp_d, fp_c, fp_b, fp_a, fp_r 170 : write(*,*) 'State in x/t' 171 : write(*,*) 'Rgas_z, rho_z, T_z, p_z, u_z' 172 : write(*,*) Rgas_z, rho_z, T_z, p_z, u_z 173 : endif 174 : C 175 : C**** Computation of thermal energy 176 : C 177 : do iord = 0, nordpo , 1 178 : acv_z(iord) = ((fp_l + fp_d) * acv_l(iord)) + 179 : & ((fp_c + fp_b) * acv_b(iord)) + 180 : & ((fp_a + fp_r) * acv_r(iord)) 181 : enddo 182 : C 182 : call prith1(nordpo, acv_z, Tmaxcv, T_z, e_z, cv_z) 183 : C 184 : hf_z = ((fp_l + fp_d) * hf_l) + ((fp_c + fp_b) * hf_b) 185 : & + ((fp_a + fp_r) * hf_r) 186 : ecin_z = 0.5d0 * u_z * u_z 187 : C 188 : C**** Interfacial flux 189 : C According to NKONGA, Comput Methods Appl. Mech Engnr 190, 2000 190 : C \dep{U}{x} + \dep{F(U)}{x} = 0 191 : C z = speed of the surface 192 : C f_z = F - z U 193 : C 194 : flurho = rho_z * (u_z - z) 195 : C write(*,*) 'flurho', flurho 196 : fluru = (flurho * u_z) + p_z 197 : fluret = (flurho * (e_z + hf_z + 198 : & ecin_z)) + (p_z * u_z) 199 : C write(*,*) (et_z + ef_z + 200 : C & ecin_z), flurho 201 : C write(*,*) p_z, u_z 202 : return 203 : end 204 :
© Cast3M 2003 - All rights reserved.
Disclaimer