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