1 : C FLUFUN SOURCE BECC 09/12/07 21:15:24 6579
2 : subroutine flufun(nordpo,
3 : & Tmaxcv, acv_z,
4 : & z,
5 : & vit_l, vit_ls, u_ls, D_resh, vit_rs, vit_r,
6 : & Rgas_l, acv_l, rho_l, u_l, T_l, hf_l,
7 : & rho_ls, T_ls,
8 : & Rgas_t, acv_t, rho_t, u_t, T_t, hf_t,
9 : & rho_rs, u_rs, T_rs,
10 : & Rgas_r, acv_r, rho_r, u_r, T_r, hf_r,
11 : & flurho, fluru, fluret)
12 : C
13 : c************************************************************************
14 : c
15 : c project :
16 : c
17 : c name : flufun
18 : c
19 : c description : Reactive Euler equations for a mixture of reactive
20 : C thermally perfect gases.
21 : c
22 : c Given z = x/t, we compute the flux
23 : C F(x/t) - ((x/t) * U)
24 : C We compute density, temperature and velocity by
25 : C averaging the intermediate state
26 : C We compute the other variables using these
27 : C quantities and the EOS
28 : c
29 : c language : fortran 77
30 : c
31 : c author : a. beccantini den/dm2s/sfme/ltmf
32 : c
33 : c************************************************************************
34 : c
34 : C called by : fcolte
35 : C
36 : c
37 : c************************************************************************
38 : c
39 : c**** input
40 : c
41 : c nordp0 = polynomial order of the cv
42 : c
43 : C Tmaxcv = Tmax for cv expansion
44 : C
45 : C z = x / t
46 : C
47 : C vit_l, vit_ls, u_ls, D_resh, vit_rs, vit_r =
48 : C = wave speeds
49 : C
50 : C Rgas_l, acv_l, rho_l, u_l, T_l, hf_l =
51 : C = some quantities in l
52 : C
53 : C rho_ls, T_ls = some quantities in ls
54 : C
55 : C Rgas_t, acv_t, rho_t, u_t, T_t, hf_t =
56 : C = some quantities in rss
57 : C
58 : C rho_rs, u_rs, T_rs =
59 : C = some quantities in rss
60 : C
61 : C Rgas_r, acv_r, rho_r, u_r, T_r, hf_r =
62 : C = some quantities in r
63 : C
64 : C acv_z = work vector
65 : C
66 : C**** output
67 : C
68 : c flurho, flurhou, fluret = fluxes for rho, rho u, rho e_t
69 : C
70 : c
71 : c************************************************************************
72 : c
73 : C 07/12/2009 created
74 : C
75 : c************************************************************************
76 : c
77 : c n.b.: all variables are declared
78 : c
79 : C implicit none
80 : integer nordpo, iord
81 : real*8 Tmaxcv, acv_z(0:nordpo), z
82 : & , vit_l, vit_ls, u_ls, D_resh, vit_rs, vit_r
83 : & , Rgas_l, acv_l(0:nordpo), rho_l, u_l, T_l, hf_l
84 : & , rho_ls, T_ls
85 : & , rho_t, u_t, T_t, hf_t
86 : & , Rgas_t, acv_t(0:nordpo), rho_rs, u_rs, T_rs
87 : & , Rgas_r, acv_r(0:nordpo), rho_r, u_r, T_r, hf_r
88 : & , flurho, fluru, fluret
89 : real*8 fp_l, fp_ls, fp_t, fp_rs, fp_r
90 : real*8 Rgas_z, rho_z, T_z, p_z, u_z, et_z
91 : & , ef_z, ecin_z
92 : real*8 T1, funt, cvtot, dcv
93 : logical logdeb
94 : parameter (logdeb = .false.)
95 : C
96 : C**** Weight functions
97 : C
98 : fp_l = 0.0d0
99 : fp_ls = 0.0d0
100 : fp_t = 0.0d0
101 : fp_rs = 0.0d0
102 : fp_r = 0.0d0
103 : C
104 : C This can be simplified from a computational point of view
105 : C
106 : if (z .lt. vit_l) then
107 : fp_l = 1.0d0
108 : else
109 : if(z .lt. vit_ls) then
110 : fp_l = (vit_ls - z) / (vit_ls - vit_l)
111 : fp_ls = 1.0d0 - fp_l
112 : else
113 : if(z .lt. u_ls)then
114 : fp_ls = 1.0d0
115 : else
116 : if(z .lt. D_resh)then
117 : fp_t = 1.0d0
118 : else
119 : if(z .lt. vit_rs)then
120 : fp_rs = 1.0d0
121 : else
122 : if(z .lt. vit_r)then
123 : fp_r = (z - vit_rs) / (vit_r - vit_rs)
124 : fp_rs = 1.0d0 - fp_r
125 : else
126 : fp_r = 1.0d0
127 : endif
128 : endif
129 : endif
130 : endif
131 : endif
132 : endif
133 : C
134 : C**** Shock-shock + (e.f) intermediate solution
135 : C
136 : Rgas_z = ((fp_l + fp_ls) * Rgas_l) +
137 : & (fp_t * Rgas_t) +
138 : & ((fp_rs + fp_r) * Rgas_r)
139 : rho_z = (fp_l * rho_l) + (fp_ls * rho_ls) +
140 : & (fp_t * rho_t) + (fp_rs * rho_rs) +
141 : & (fp_r * rho_r)
142 : T_z = (fp_l * T_l) + (fp_ls * T_ls) +
143 : & (fp_t * T_t) + (fp_rs * T_rs) +
144 : & (fp_r * T_r)
145 : p_z = Rgas_z * rho_z * T_z
146 : u_z = (fp_l * u_l) + (fp_ls * u_ls) +
147 : & (fp_t * u_t) + (fp_rs * u_rs) +
148 : & (fp_r * u_r)
149 : if (logdeb) then
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_ls) * acv_l(iord)) +
159 : & (fp_t * acv_t(iord)) +
160 : & ((fp_rs + fp_r) * acv_r(iord))
161 : enddo
162 : C T1 = min (T, Tmaxcv)
163 : T1 = 0.5d0 * ((T_z + Tmaxcv) - abs (T_z - Tmaxcv))
164 : C
165 : funt = 1.0D0
166 : cvtot = acv_z(0)
167 : et_z = cvtot * T1
168 : do iord = 1, nordpo, 1
169 : funt = funt * T1
170 : dcv = acv_z(iord) * funt
171 : cvtot = cvtot + dcv
172 : et_z = et_z + (dcv * T1 / (iord + 1))
173 : enddo
174 : C
175 : C if ( T > Tmaxcv ) => T1 = (T - Tmaxcv)
176 : C else T1 = 0
177 : C i.e.
178 : C T1 = max ((T - Tmaxcv),0)
179 : C
180 : T1 = 0.5d0 * ((T_z - Tmaxcv) + abs (T_z - Tmaxcv))
181 : et_z = et_z + (cvtot * T1)
182 : C
183 : ef_z = ((fp_l + fp_ls) * hf_l) + (fp_t * hf_t)
184 : & + ((fp_r + fp_rs) * hf_r)
185 : ecin_z = 0.5d0 * u_z * u_z
186 : C write(*,*) 'et, ecin, ef', et_z, ecin_z, ef_z
187 : C write(*,*) 'rho_z, u_z, p_z', rho_z, u_z, p_z
188 : C
189 : C**** Interfacial flux
190 : C According to NKONGA, Comput Methods Appl. Mech Engnr 190, 2000
191 : C \dep{U}{x} + \dep{F(U)}{x} = 0
192 : C z = speed of the surface
193 : C f_z = F - z U
194 : C
195 : flurho = rho_z * (u_z - z)
196 : C write(*,*) 'flurho', flurho
197 : fluru = (flurho * u_z) + p_z
198 : fluret = (flurho * (et_z + ef_z +
199 : & ecin_z)) + (p_z * u_z)
200 : C write(*,*) (et_z + ef_z +
201 : C & ecin_z), flurho
202 : C write(*,*) p_z, u_z
203 : return
204 : end
205 :
© Cast3M 2003 - All rights reserved.
Disclaimer