C FUNRE1 SOURCE BECC 09/12/07 21:15:31 6579 subroutine funre1(nordpo, tmaxcv, rgas_a, acv_a, & rgas_b, acv_b, q, k0, T_ag, T_bg, & r_a, p_a, u_a, T_a, & r_b, p_b, u_b, T_b, Dsw, logan) c************************************************************************ c c projet : c c name : funre1 c c description : computation of the weak state behind a reactive c shock c c language : fortran 77 c c author : a. beccantini den/dm2s/sfme/ltmf c c************************************************************************ c c callen by : racre2 c c called program : none c c************************************************************************ c c c**** input c c nordpo = polynomials degree c tmaxcv = threshold temperature for cv(t) c rgas_a = r of the gas in r and rs c acv_a = \sum_{i=1,nesp} y_i*acv_a{i,j} c rgas_b = r of the gas in rss c acv_b = \sum_{i=1,nesp} y_i*acv_b{i,j} c q = realeased heat per mass unit c k0 = fundamental flame speed C T_ag,T_bg = guess temperatures for a and b C (to compute properties) C r_a = density in rs C p_a = pressure in rs c u_a = velocity in rs C T_a = temperature in rs c**** output c r_b = approximate density in rss (behind the shock) c p_b = approximate pressure in rss c u_b = approximate speed in rss C T_b = approximate temperature in rss C (NB if T_bg = T_b and T_ag = T_a, rss thus C computed is exact) C Dsw = approximate shock wave speed... it could C differs from u_a + k0 if (k0lim < k0) C Nevertheless, if logsto = .true., if k0lim C < k0 we stop! c logan = an anomaly is detected c c************************************************************************ c c created the 07/12/2009 c c implicit none integer nordpo real*8 tmaxcv, rgas_a, acv_a(1:(nordpo+1)) & ,rgas_b, acv_b(1:(nordpo+1)) real*8 q, k0 real*8 T_ag, T_bg real*8 r_a, p_a, u_a, T_a real*8 r_b, p_b, u_b, T_b, Dsw real*8 gam_a, e_a, cv_a, gam_b, e_b, cv_b real*8 ttq, k02ad, h_a real*8 a, b, c, delta, delta1, k0cj2, k0new real*8 csi logical logan CC CC**** We write the input CC C write(*,*) 'We enter funre1' C C**** Computation of gam_a and gam_b C call prith1(nordpo, acv_a, Tmaxcv, T_ag, e_a, cv_a) gam_a = (cv_a + Rgas_a) / cv_a call prith1(nordpo, acv_b, Tmaxcv, T_bg, e_b, cv_b) gam_b = (cv_b + Rgas_b) / cv_b ttq = q + ((cv_b * T_bg) - e_b) - & ((cv_a * T_ag) - e_a) C h_a and h_b computed as we deal with cp gases h_a = (gam_a / (gam_a - 1.0D0)) * Rgas_a * T_a C r_a = P_a / (Rgas_a * T_a) C C**** Computation of r_b C k02ad = k0 * k0 / (Rgas_a * T_a) if (k02ad .lt. 1.0D0)then C if ( .true. ) then C if ( .false. ) then C write(*,*) 'funre1.f. Parte 1' c = k02ad b = k02ad + 1.0D0 b = b * gam_b / (gam_b + 1.0D0) a = 2.0D0 * (h_a + ttq) / (Rgas_a * T_a) a = a + k02ad a = a * (gam_b - 1.0D0) / (gam_b + 1.0D0) delta = (b * b) - (a * c) if (delta .lt. 0.0D0) then C if ( .true. ) then C write(*,*) 'funre1.f. Parte 1.2' C We are in CJDF regime a = 1.0D0 / (gam_b * gam_b) b = 1.0D0 - a b = b * (h_a + ttq) b = b - (Rgas_a * T_a) c = (Rgas_a * T_a)**2 delta1 = (b * b) - (a * c) if (delta1 .lt. 0.0D0) then write(*,*) 'subroutine funre1.f' write(*,*) 'We cannot compute k0cj' write(*,*) 'ANOMALY DETECTED' logan = .true. goto 9999 endif k0cj2 = (b - sqrt(delta1)) / a if (k0cj2 .lt. 0.0D0) then write(*,*) 'subroutine funre1.f' write(*,*) 'k0cj^2 < 0' write(*,*) 'ANOMALY DETECTED' logan = .true. goto 9999 endif k02ad = k0cj2 / (Rgas_a * T_a) b = k02ad + 1.0D0 b = b * gam_b / (gam_b + 1.0D0) a = 2.0D0 * (h_a + ttq) / (Rgas_a * T_a) a = a + k02ad a = a * (gam_b - 1.0D0) / (gam_b + 1.0D0) csi = b / a k0new = k0cj2 ** 0.5D0 else k0new = k0 csi = b + sqrt(delta) csi = csi / a endif r_b = r_a * csi else C write(*,*) 'funre1.f. Parte 2' c = 2.0D0 * (h_a + ttq) c = c / (k0 * k0) c = c + 1.0D0 c = c * (gam_b - 1.0D0) / (gam_b + 1.0D0) b = Rgas_a * T_a / (k0 * k0) b = b + 1.0D0 b = b * gam_b / (gam_b + 1.0D0) delta = (b * b) - c if (delta .lt. 0.0D0) then C if ( .true. ) then C write(*,*) 'funre1.f. Parte 2.2' a = 1.0D0 / (gam_b * gam_b) b = 1.0D0 - a b = b * (h_a + ttq) b = b - (Rgas_a * T_a) c = (Rgas_a * T_a) c = c * c delta1 = (b * b) - (a * c) if (delta1 .lt. 0.0D0) then write(*,*) 'subroutine funre1.f' write(*,*) 'We cannot compute k0cj' write(*,*) 'ANOMALY DETECTED' logan = .true. goto 9999 endif k0cj2 = (b - sqrt(delta1)) / a if (k0cj2 .lt. 0.0D0) then write(*,*) 'subroutine funre1.f' write(*,*) 'k0cj^2 < 0' write(*,*) 'ANOMALY DETECTED' logan = .true. goto 9999 endif c = 2.0D0 * (h_a + ttq) c = c / k0cj2 c = c + 1.0D0 c = c * (gam_b - 1.0D0) / (gam_b + 1.0D0) b = Rgas_a * T_a / k0cj2 b = b + 1.0D0 b = b * gam_b / (gam_b + 1.0D0) delta1 = b * b - c C if (abs (delta1) > 1.0D-6) then if (abs (delta1) .gt. 1.0D-6) then write(*,*) 'subroutine funre1.f' write(*,*) 'CJDF' write(*,*) 'delta != 0' write(*,*) 'ANOMALY DETECTED' logan = .true. goto 9999 endif csi = b k0new = k0cj2 ** 0.5D0 else csi = b - sqrt (delta) k0new = k0 endif r_b = r_a / csi endif csi = r_a / r_b Dsw = k0new + u_a u_b = Dsw - (k0new * csi) p_b = p_a + (r_a * k0new * k0new * (1.0d0 - csi)) if (p_b .lt. 0.0D0) then write(*,*) 'subroutine funre1.f' write(*,*) 'p_b < 0' write(*,*) 'ANOMALY DETECTED' logan = .true. goto 9999 endif T_b = P_b / (Rgas_b * r_b) C write(*,*) 'r_b, p_b, u_b, T_b, Dsw, k0new' C write(*,'(6E12.4)') r_b, p_b, u_b, T_b, Dsw C write(*,*) r_b, p_b, u_b, T_b, Dsw, k0new C stop CC write(*,*) 'We exit funre1' CC 9999 continue return end