OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps41.F
Go to the documentation of this file.
1Copyright> OpenRadioss
2Copyright> Copyright (C) 1986-2025 Altair Engineering Inc.
3Copyright>
4Copyright> This program is free software: you can redistribute it and/or modify
5Copyright> it under the terms of the GNU Affero General Public License as published by
6Copyright> the Free Software Foundation, either version 3 of the License, or
7Copyright> (at your option) any later version.
8Copyright>
9Copyright> This program is distributed in the hope that it will be useful,
10Copyright> but WITHOUT ANY WARRANTY; without even the implied warranty of
11Copyright> MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12Copyright> GNU Affero General Public License for more details.
13Copyright>
14Copyright> You should have received a copy of the GNU Affero General Public License
15Copyright> along with this program. If not, see <https://www.gnu.org/licenses/>.
16Copyright>
17Copyright>
18Copyright> Commercial Alternative: Altair Radioss Software
19Copyright>
20Copyright> As an alternative to this open-source version, Altair also offers Altair Radioss
21Copyright> software under a commercial license. Contact Altair to discuss further if the
22Copyright> commercial version may interest you: https://www.altair.com/radioss/.
23!||====================================================================
24!|| sigeps41 ../engine/source/materials/mat/mat041/sigeps41.F
25!||--- called by ------------------------------------------------------
26!|| mulaw ../engine/source/materials/mat_share/mulaw.F90
27!|| mulaw8 ../engine/source/materials/mat_share/mulaw8.F90
28!||--- calls -----------------------------------------------------
29!|| burn_fraction_ireac1 ../engine/source/materials/mat/mat041/sigeps41.F
30!|| burn_fraction_ireac2 ../engine/source/materials/mat/mat041/sigeps41.F
31!|| mixture_equilibrium ../engine/source/materials/mat/mat041/sigeps41.F
32!||--- uses -----------------------------------------------------
33!|| output_mod ../common_source/modules/output/output_mod.F90
34!|| root_finding_algo_mod ../common_source/modules/root_finding_algo_mod.F90
35!||====================================================================
36 SUBROUTINE sigeps41(OUTPUT,
37 1 NEL ,NUPARAM,NUVAR ,
38 2 TIME ,TIMESTEP,UPARAM ,RHO0 ,RHO ,
39 3 VOLUME ,EINT ,
40 7 SIGOXX ,SIGOYY ,SIGOZZ ,
41 8 SIGNXX ,SIGNYY ,SIGNZZ ,SIGNXY ,SIGNYZ ,SIGNZX ,
42 9 SIGVXX ,SIGVYY ,SIGVZZ ,SIGVXY ,SIGVYZ ,SIGVZX ,
43 A SOUNDSP,VISCMAX,UVAR ,DELTVOL ,BFRAC ,TEMP )
44C-----------------------------------------------
45C D e s c r i p t i o n
46C-----------------------------------------------
47 USE root_finding_algo_mod
48 USE output_mod, ONLY : output_
49C-----------------------------------------------
50C D e s c r i p t i o n
51C-----------------------------------------------
52C MULAW
53C |__SIGEPS41 (LEE TARVER)
54C |__ MIX (JWL EOS solve mixture state for both reagent & product, Pr=Pp=P).
55C | iterative solver
56C | |__ JWL_EOS_STATE (JWL EoS single phase / reagent)
57C | |__ JWL_EOS_STATE (JWL EoS single phase / product)
58C | |__ FTNCH
59C |
60C |__ BURN_FRACTION_IREAC1 (update F function, case Ireac=1)
61C |__ BURN_FRACTION_IREAC2 (update F function, case Ireac=2)
62C
63C ________________________________
64C IREAC FLAG
65C
66C *** IREAC = 1 ***
67C ORIGINAL 2-TERM-MODEL --> dF/dT = IGNITION TERM + GROWTH TERM
68C IGNITION TERM : ratei = I * (1.-F)**b * (etar-1.)**x , if P>0
69C GROWTH TERM : rateg = G1 * (1.-F)**b * F**d * pres**y , if P>0
70C
71C *** IREAC = 2 ***
72C EXTENDED 3-TERM-MODEL --> dF/dT = IGNITION TERM + GROWTH TERM 1 + GROWTH TERM 2
73C IGNITION TERM : ratei = I * (1.-F+tol)**b * abs(etar-1.-ccrit)**x , if 0<f<Figmax & P>0 & compression threshold (etar-1>=ccrit)
74C GROWTH TERM 1 : rateg1 = G1 * (1.-F+tol)**c * (F+tol)**d * pres**y , if 0<f<FG1max & P>0
75C GROWTH TERM 2 : rateg2 = G2 * (1.-F+tol)**e * (F+tol)**g * pres**z , if FG2min<f<1 & P>0
76C
77C ________________________________
78C NOTATION
79C index 'p' : for 'p'roducts of detonation
80C index 'r' : for 'r'eagent (unreacted explosive)
81C fc : Function F described by Lee-Tarver model dF/dt = ...
82C
83C-----------------------------------------------
84C I m p l i c i t T y p e s
85C-----------------------------------------------
86#include "implicit_f.inc"
87#include "comlock.inc"
88C-----------------------------------------------
89C G l o b a l P a r a m e t e r s
90C-----------------------------------------------
91#include "com01_c.inc"
92#include "com06_c.inc"
93C---------+---------+---+---+--------------------------------------------
94C VAR | SIZE |TYP| RW| DEFINITION
95C---------+---------+---+---+--------------------------------------------
96C NEL | 1 | I | R | SIZE OF THE ELEMENT GROUP NEL
97C NUPARAM | 1 | I | R | SIZE OF THE USER PARAMETER ARRAY
98C NUVAR | 1 | I | R | NUMBER OF USER ELEMENT VARIABLES
99C---------+---------+---+---+--------------------------------------------
100C---------+---------+---+---+--------------------------------------------
101C TIME | 1 | F | R | CURRENT TIME
102C TIMESTEP| 1 | F | R | CURRENT TIME STEP
103C UPARAM | NUPARAM | F | R | USER MATERIAL PARAMETER ARRAY
104C RHO0 | NEL | F | R | INITIAL DENSITY
105C RHO | NEL | F | R | DENSITY
106C VOLUME | NEL | F | R | VOLUME
107C EINT | NEL | F | R | TOTAL INTERNAL ENERGY
108C EPSPXX | NEL | F | R | STRAIN RATE XX
109C EPSPYY | NEL | F | R | STRAIN RATE YY
110C ... | | | |
111C EPSXX | NEL | F | R | STRAIN XX
112C EPSYY | NEL | F | R | STRAIN YY
113C ... | | | |
114C SIGOXX | NEL | F | R | OLD ELASTO PLASTIC STRESS XX
115C SIGOYY | NEL | F | R | OLD ELASTO PLASTIC STRESS YY
116C ... | | | |
117C---------+---------+---+---+--------------------------------------------
118C SIGNXX | NEL | F | W | NEW ELASTO PLASTIC STRESS XX
119C SIGNYY | NEL | F | W | NEW ELASTO PLASTIC STRESS YY
120C ... | | | |
121C SIGVXX | NEL | F | W | VISCOUS STRESS XX
122C SIGVYY | NEL | F | W | VISCOUS STRESS YY
123C ... | | | |
124C SOUNDSP | NEL | F | W | SOUND SPEED (NEEDED FOR TIME STEP)
125C VISCMAX | NEL | F | W | MAXIMUM DAMPING MODULUS(NEEDED FOR TIME STEP)
126C---------+---------+---+---+--------------------------------------------
127C UVAR |NEL*NUVAR| F |R/W| USER ELEMENT VARIABLE ARRAY
128C---------+---------+---+---+--------------------------------------------
129C DELTVOL | NEL | F |R | VOLUME VARIATION
130C---------+---------+---+---+--------------------------------------------
131C-----------------------------------------------
132C D u m m y A r g u m e n t s
133C-----------------------------------------------
134 TYPE(output_), INTENT(INOUT) :: OUTPUT
135 INTEGER NEL, NUPARAM, NUVAR
136
137 my_real, INTENT(IN) ::
138 . TIME,TIMESTEP,UPARAM(NUPARAM),
139 . RHO(NEL),RHO0(NEL),VOLUME(NEL),
140 . SIGOXX(NEL),SIGOYY(NEL),SIGOZZ(NEL),
141 . deltvol(nel)
142
143 my_real, INTENT(INOUT) ::
144 . signxx(nel),signyy(nel),signzz(nel),
145 . signxy(nel),signyz(nel),signzx(nel),
146 . sigvxx(nel),sigvyy(nel),sigvzz(nel),
147 . sigvxy(nel),sigvyz(nel),sigvzx(nel),
148 . soundsp(nel),viscmax(nel),eint(nel),bfrac(nel),temp(nel)
149
150 my_real,INTENT(INOUT) :: uvar(nel,nuvar)
151C-----------------------------------------------
152C L o c a l V a r i a b l e s
153C-----------------------------------------------
154 my_real ar,br,r1r,r2r,r3r,cvr,etar,
155 . ap,bp,r1p,r2p,r3p,cvp,etap,
156 . epsil,ftol, tmp,
157 . i_,b_,x_,g1,d_,y_,g2,
158 . z_,g_,c_,e_,figmax,fg1max,fg2min,
159 . cappa,chi,ccrit,tol,cv,dedv,
160 . artvisc,dt,pres,fc,cburn,enq,shear,vol_rel,dvol_rel
161 my_real oldfc,tem1,tem2,chydro,chydro2,dpdmu,dpde,heat
162 my_real mt,pold,er,wfextt,wr,wp_coef
163 integer itrmax,i_reac,i
164 my_real :: beta !volumic fraction of reagent
165 my_real, dimension(25) :: funct_parameter
166
167C-----------------------------------------------
168C S o u r c e L i n e s
169C-----------------------------------------------
170 wfextt=zero
171 if (time == zero) then
172 do i=1,nel
173 uvar(i,1) = zero ! Burn Fraction (F)
174 uvar(i,2) = uparam(39) ! temperature
175 uvar(i,3) = uparam(14) ! cv mixture, volumetric heat capacity (SI: J/m3/K)
176 uvar(i,4) = one ! relative volume
177 uvar(i,5) = zero ! dedv
178 uvar(i,6) = zero ! eta_p (eta:=rho/rho0)
179 uvar(i,7) = one ! eta_r (eta:=rho/rho0)
180 enddo
181 endif
182
183 i_reac = int(uparam(1))
184 itrmax = int(uparam(18))
185 dt = timestep
186 ar = uparam(2)
187 br = uparam(3)
188 r1r = uparam(4)
189 r2r = uparam(5)
190 r3r = uparam(6)
191 wr = uparam(7)
192 cvr = uparam(14)
193 ap = uparam(8)
194 bp = uparam(9)
195 r1p = uparam(10)
196 r2p = uparam(11)
197 r3p = uparam(12)
198 wp_coef = uparam(13)
199 cvp = uparam(15)
200 enq = uparam(16)
201 epsil = uparam(17)
202 ftol = uparam(19)
203 i_ = uparam(20)
204 b_ = uparam(21)
205 x_ = uparam(22)
206 g1 = uparam(23)
207 d_ = uparam(24)
208 y_ = uparam(25)
209 cappa = uparam(26)
210 chi = uparam(27)
211 tol = uparam(28)
212 ccrit = uparam(29)
213 g2 = uparam(30)
214 c_ = uparam(31)
215 e_ = uparam(32)
216 g_ = uparam(33)
217 z_ = uparam(34)
218 figmax = uparam(35)
219 fg1max = uparam(36)
220 fg2min = uparam(37)
221 shear = uparam(38)
222c
223
224 funct_parameter(1:25) = zero
225 funct_parameter(1) = ar
226 funct_parameter(2) = br
227 funct_parameter(3) = r1r
228 funct_parameter(4) = r2r
229 funct_parameter(5) = r3r
230 funct_parameter(6) = cvr
231 ! funct_parameter(7) = tmp output
232 ! funct_parameter(8) = dedvr output
233 ! funct_parameter(9) = preac output
234 ! funct_parameter(10) = bthr ! output
235 ! funct_parameter(11) = dpdtr! output
236 ! funct_parameter(12) = enr! output
237
238 funct_parameter(13) = ap
239 funct_parameter(14) = bp
240 funct_parameter(15) = r1p
241 funct_parameter(16) = r2p
242 funct_parameter(17) = r3p
243 funct_parameter(18) = cvp
244 ! funct_parameter(19) = dedvp
245 ! funct_parameter(20) = pprod
246 ! funct_parameter(21) = bthp
247 ! funct_parameter(22) = dpdtp
248 ! funct_parameter(23) = enp
249 ! funct_parameter(24) = etac
250 ! funct_parameter(25) = fc
251
252
253 do i=1,nel
254 vol_rel = rho0(i) / rho(i) !no unit
255 pres = -(sigoxx(i)+sigoyy(i)+sigozz(i)) * third ! pressure
256 pold = pres
257 artvisc = zero !pseudo viscosity
258 mt=rho(i)*volume(i) !mass
259 IF (time > zero .AND. iale + ieuler > 0) THEN
260 uvar(i, 1) = uvar(i, 1) / mt
261 uvar(i, 1) = max(zero, min(uvar(i, 1), one))
262 ENDIF
263 fc = uvar(i,1) ! F function (no unit)
264 cv = uvar(i,3) ! volumetric heat capacity (SI: J/m3/K)
265 tmp = uparam(39) ! temperature K
266 if(time > zero .AND. mt > em20) tmp = eint(i)/(mt*cv)
267 dvol_rel = vol_rel - uvar(i,4)
268 dedv = uvar(i,5)
269 etap = uvar(i,6)
270 etar = uvar(i,7)
271C=======================================================================
272c Pressure and Sound Speed
273c
274c initi
275 oldfc = fc
276c
277c hydrodynamic state (pass 1)
278 tem1 = (pres+dedv+artvisc) / cv ! (Cv SI:J/m3/K = Pa/K)
279 tmp = tmp - tem1*dvol_rel
280 heat=zero
281 dpdmu=zero
282 dpde=zero
283 CALL mixture_equilibrium(ar ,br ,r1r ,r2r ,r3r ,cvr ,etar,
284 . ap ,bp ,r1p ,r2p ,r3p ,cvp ,etap,
285 . dpdmu,dpde,ftol,enq,
286 . tmp,heat,vol_rel ,fc,cv,dedv,pres, beta,funct_parameter)
287c
288c hydrodynamic state (pass 2)
289 tem2 = (pres+dedv+artvisc) / cv ! (Cv SI:J/m3/K = Pa/K)
290 tmp = tmp - half*(tem2-tem1)*dvol_rel
291 CALL mixture_equilibrium(ar ,br ,r1r ,r2r ,r3r ,cvr ,etar,
292 . ap ,bp ,r1p ,r2p ,r3p ,cvp ,etap,
293 . dpdmu,dpde,ftol,enq,
294 . tmp,heat,vol_rel ,fc,cv,dedv,pres, beta,funct_parameter)
295c
296c dF/dt (F function)
297 !ORIGINAL 2-TERM-MODEL
298 if (i_reac == 1) CALL burn_fraction_ireac1(i_,b_,x_,g1,d_,y_,cappa,artvisc,dt,pres,etar,fc,cburn,oldfc)
299 !EXTENDED 3-TERM-MODEL
300 if (i_reac == 2) CALL burn_fraction_ireac2(i_,b_,x_,g1,d_,y_,g2,z_,g_,c_,e_,figmax,fg1max,fg2min,
301 . cappa,chi,ccrit,tol,artvisc,dt,pres,etar,fc,cburn,oldfc)
302 tmp = tmp + (fc - oldfc) * heat/cv
303c sound speed (time step criteria)
304 chydro2 = dpdmu + abs(pres)*(vol_rel*vol_rel)*dpde + onep33*shear
305 chydro2=max(chydro2,max(wr+one,wp_coef+one)*abs(pres)/max(em30,rho(i)))
306 chydro = sqrt(chydro2)
307c
308C BUFFER STORAGE
309C --------------
310 soundsp(i) = max(chydro,cburn)
311 viscmax(i) = zero
312c stress tensor
313 signxx(i) = -pres
314 signyy(i) = -pres
315 signzz(i) = -pres
316 signxy(i) = zero
317 signyz(i) = zero
318 signzx(i) = zero
319c viscous stress tensor
320 sigvxx(i) = zero
321 sigvyy(i) = zero
322 sigvzz(i) = zero
323 sigvxy(i) = zero
324 sigvyz(i) = zero
325 sigvzx(i) = zero
326c material buffer
327 uvar(i,1) = fc
328 uvar(i,2) = tmp
329 uvar(i,3) = cv
330 uvar(i,4) = vol_rel
331 uvar(i,5) = dedv
332 uvar(i,6) = etap !products:rho_p/rho0_p
333 uvar(i,7) = etar !reagents:rho_r/rho0_r
334C Burn fraction
335 bfrac(i) = fc ! burn fraction (function F)
336 uvar(i,8) = one - beta !volume fraction of products
337C energie released by the reaction = D(energie totale) - (-P.DV)
338 er=mt*cv*tmp + half*(pres+pold)*deltvol(i) !misunderstanding between volumetric heat capacity (J/m3/K = Pa/K) and heat capacity(J/K) ? law41 to be verified.
339 wfextt=wfextt+er-eint(i)
340 eint(i)=mt*cv*tmp + half*(pres+pold)*deltvol(i)
341C EINT(I)=EINT(I)-0.5*(pres+pold)*DELTVOL(I) done in mulaw.
342 temp(i) = tmp !temperature
343 enddo
344c
345!$OMP ATOMIC
346 output%TH%WFEXT=output%TH%WFEXT+wfextt
347C-----------------------------------------------
348 RETURN
349 END
350C-----------------------------------------------
351
352C-----------------------------------------------
353!||====================================================================
354!|| burn_fraction_ireac1 ../engine/source/materials/mat/mat041/sigeps41.F
355!||--- called by ------------------------------------------------------
356!|| sigeps41 ../engine/source/materials/mat/mat041/sigeps41.F
357!||====================================================================
358 SUBROUTINE burn_fraction_ireac1(I_,B_,X_,G1,D_,Y_,cappa,artvisc,dt,pres,etar,fc,cburn,oldfc)
359C-----------------------------------------------
360C D e s c r i p t i o n
361C-----------------------------------------------
362
363C *** IREAC = 1 ***
364C ORIGINAL 2-TERM-MODEL --> dF/dT = IGNITION TERM + GROWTH TERM
365C IGNITION TERM : ratei = I * (1.-F)**b * (etar-1.)**x , if P>0
366C GROWTH TERM : rateg = G1 * (1.-F)**b * F**d * pres**y , if P>0
367c
368c bounds :
369c -dfc limited during 1 cycle (chi)
370c -dfc limited within a shock front (cappa)
371c cburn not calculated here
372c
373c input : etar, pres, oldfc, dt
374c output : fc, cburn
375c
376C-----------------------------------------------
377C I m p l i c i t T y p e s
378C-----------------------------------------------
379#include "implicit_f.inc"
380C-----------------------------------------------
381C D u m m y A r g u m e n t s
382C-----------------------------------------------
383 my_real
384 . i_,b_,x_,g1,d_,y_,cappa,
385 . artvisc,dt,pres,etar,fc,cburn,oldfc
386C-----------------------------------------------
387C L o c a l V a r i a b l e s
388C-----------------------------------------------
389 my_real ratei,rateg,dfci,dfcg
390C-----------------------------------------------
391C P r e c o n d i t i o n s
392C-----------------------------------------------
393 cburn = zero
394 if (pres <= zero) return !no rate to calculate for material expansion
395 if (artvisc >= (cappa*pres)) return !
396 if (fc == one) return !everything already reacted, 100% products, Fc remains 1.00
397C-----------------------------------------------
398C S o u r c e L i n e s
399C-----------------------------------------------
400c ignition term
401 ratei = i_ * (one -fc)**b_ * (etar-one)**x_
402 if ((etar - one) <= zero) ratei=zero
403 dfci = ratei * dt
404c growth term
405 rateg = g1 * (one - fc)**b_ * fc**d_ * pres**y_
406 dfcg = rateg * dt
407c F function (0. < F < 1.)
408 fc = oldfc + max(dfci,zero) + max(dfcg,zero)
409 if (fc >= one) fc = one
410C-----------------------------------------------
411 RETURN
412 END SUBROUTINE burn_fraction_ireac1
413C-----------------------------------------------
414
415
416!||====================================================================
417!|| burn_fraction_ireac2 ../engine/source/materials/mat/mat041/sigeps41.F
418!||--- called by ------------------------------------------------------
419!|| sigeps41 ../engine/source/materials/mat/mat041/sigeps41.F
420!||====================================================================
421 SUBROUTINE burn_fraction_ireac2(I_,B_,X_,G1,D_,Y_,G2,Z_,G_,C_,E_,Figmax,FG1max,FG2min,cappa,chi,ccrit,tol,
422 . artvisc,dt,pres,etar,fc,cburn,oldfc)
423C-----------------------------------------------
424C D e s c r i p t i o n
425C-----------------------------------------------
426C *** IREAC = 2 ***
427C EXTENDED 3-TERM-MODEL --> dF/dT = IGNITION TERM + GROWTH TERM 1 + GROWTH TERM 2
428C IGNITION TERM : ratei = I * (1.-F+tol)**b * abs(etar-1.-ccrit)**x , if 0<f<Figmax & P>0 & compression threshold (etar-1>=ccrit)
429C GROWTH TERM 1 : rateg1 = G1 * (1.-F+tol)**c * (F+tol)**d * pres**y , if 0<f<FG1max & P>0
430C GROWTH TERM 2 : rateg2 = G2 * (1.-F+tol)**e * (F+tol)**g * pres**z , if FG2min<f<1 & P>0
431c
432c bounds :
433c - dfc limited during 1 cycle (chi)
434c - dfc limited within a shock front (cappa)
435c buirn not calculated here
436c
437c input : etar, pres, oldfc, dt
438c output : fc, cburn
439c
440C-----------------------------------------------
441C I m p l i c i t T y p e s
442C-----------------------------------------------
443#include "implicit_f.inc"
444C-----------------------------------------------
445C D u m m y A r g u m e n t s
446C-----------------------------------------------
447 my_real i_,b_,x_,g1,d_,y_,g2,z_,g_,c_,e_,
448 . figmax,fg1max,fg2min,cappa,chi,ccrit,tol,
449 . artvisc,dt,pres,etar,fc,cburn,oldfc
450C-----------------------------------------------
451C L o c a l V a r i a b l e s
452C-----------------------------------------------
453 my_real ratei,rateg,rateg1,rateg2,dfci,dfcg
454C-----------------------------------------------
455C P r e c o n d i t i o n s
456C-----------------------------------------------
457 cburn = zero
458 if (pres <= zero) return
459 if (fc == one) return
460 if (artvisc >= (cappa*pres)) return
461 if (fc <= figmax) then
462C-----------------------------------------------
463C S o u r c e L i n e s
464C-----------------------------------------------
465 ! ignition term
466 ratei = i_ * (one-fc+tol)**b_ * abs(etar-one-ccrit)**x_
467 if ((etar- one -ccrit)<=zero) ratei = zero
468 dfci = ratei * dt
469 fc = oldfc + max(dfci,zero)
470 if (fc > figmax) fc = figmax
471 endif
472
473 ! growth term 1
474 rateg1 = g1 * (one -fc+tol)**c_ * (fc+tol)**d_ * pres**y_
475 if (fc > fg1max) rateg1 = zero
476 ! growth term 2
477 rateg2 = g2 * (one -fc+tol)**e_ * (fc+tol)**g_ * pres**z_
478 if (fc < fg2min) rateg2 = zero
479 ! F function (0. < F < 1.)
480 rateg = rateg1 + rateg2
481 dfcg = rateg * dt
482 fc = fc + max(dfcg,zero)
483 if ((fc-oldfc)>chi) fc = oldfc + chi
484 if (fc >= one) fc = one
485C-----------------------------------------------
486 RETURN
487 END SUBROUTINE burn_fraction_ireac2
488C-----------------------------------------------
489
490
491!||====================================================================
492!|| mixture_equilibrium ../engine/source/materials/mat/mat041/sigeps41.F
493!||--- called by ------------------------------------------------------
494!|| sigeps41 ../engine/source/materials/mat/mat041/sigeps41.F
495!||--- calls -----------------------------------------------------
496!|| jwl_eos_delta ../engine/source/materials/mat/mat041/jwl_eos_mod.F90
497!|| jwl_eos_state ../engine/source/materials/mat/mat041/jwl_eos_mod.F90
498!||--- uses -----------------------------------------------------
499!|| jwl_eos_mod ../engine/source/materials/mat/mat041/jwl_eos_mod.F90
500!|| root_finding_algo_mod ../common_source/modules/root_finding_algo_mod.F90
501!||====================================================================
502 SUBROUTINE mixture_equilibrium(
503 . ar ,br ,r1r ,r2r ,r3r ,cvr ,etar,
504 . ap ,bp ,r1p ,r2p ,r3p ,cvp ,etap,
505 . dpdmu,dpde ,ftol,enq ,
506 . tmp ,heat ,vol_rel,fc,cv,dedv ,pres, beta,funct_parameter)
507 USE root_finding_algo_mod
508 USE jwl_eos_mod, ONLY: jwl_eos_delta, jwl_eos_state
509C-----------------------------------------------
510C D e s c r i p t i o n
511C-----------------------------------------------
512c fc massic fraction (F function) ;
513c beta : volume fraction of reagent (unreacted explosive)
514c T,v,fc known : find beta such as Pr=Pp (hypothesis : equilibrium between two phase for system closure)
515c Iterative method : iterate on b value until delta = Pp-Pr < eps
516c
517c input : vol_rel, etar, etap, fc, tmp
518c output : etar, etap, pres, dedv,cv,heat,dpde,dpdmu
519c local : beta
520c
521C-----------------------------------------------
522C I m p l i c i t T y p e s
523C-----------------------------------------------
524#include "implicit_f.inc"
525C-----------------------------------------------
526C D u m m y A r g u m e n t s
527C-----------------------------------------------
528 my_real ar ,br ,r1r ,r2r ,r3r ,cvr ,etar,
529 . ap ,bp ,r1p ,r2p ,r3p ,cvp ,etap,
530 . ftol,heat,dpdmu,dpde,enq, tmp,
531 . vol_rel ,fc,cv,dedv,pres, beta
532 my_real, dimension(25), intent(inout) :: funct_parameter
533C-----------------------------------------------
534C L o c a l V a r i a b l e s
535C-----------------------------------------------
536 integer iter,i1
537 my_real delta,fc1,etac,guess,pmax,
538 . trans,transr,transp,trans1,trans2,dbdeta,detar,detap,
539 . dedr,dedp,bth,dbdt,detrdt,detpdt,dpdt,dbdf,detrdf,
540 . detpdf,cvr0,cvp0,slope,tl,tu,a,b,c,
541 . bthp,dedvp,preac,pprod,dpdtp,dpdtr,bthr,dedvr,
542 . enp,enr
543 my_real :: lower_bound,upper_bound
544C-----------------------------------------------
545C S o u r c e L i n e s
546C-----------------------------------------------
547c initi.
548 etac = one/vol_rel
549 fc1 = one -fc
550 if (fc < ftol) then
551c reaction has not started : F=0., beta=Vfrac(reagent)=1.
552 etar = fc1 * etac
553 CALL jwl_eos_state(ar,br,r1r,r2r,r3r,cvr,etar,tmp,dedvr,pres,bthr,dpdtr,enr)
554 beta = one
555 cv = cvr
556 dedv = dedvr
557 dpde = dpdtr / cv
558 dpdmu = bthr + dpdtr*dedv / (cv*etac**2)
559 etap = em6
560 CALL jwl_eos_state(ap,bp,r1p,r2p,r3p,cvp,etap,tmp,dedvp,pprod,bthp,dpdtp,enp)
561 heat = enq-enp+enr
562 return
563c
564 elseif (fc > (one -ftol)) then
565c reaction has finished : F=1., beta=Vfrac(reagent)=0.
566 etap = fc * etac
567 CALL jwl_eos_state(ap,bp,r1p,r2p,r3p,cvp,etap,tmp,dedvp,pres,bthp,dpdtp,enp )
568 cv = cvp
569 beta = zero
570 heat = enq
571 dedv = dedvp
572 dpde = dpdtp / cv
573 dpdmu = bthp + dpdtp*dedv / (cv*etac**2)
574 return
575c
576 endif
577
578
579 funct_parameter(7) = tmp
580 funct_parameter(24) = etac
581 funct_parameter(25) = fc
582
583 funct_parameter(8) = zero !dedvr output
584 funct_parameter(9) = zero !preac output
585 funct_parameter(10) = zero !bthr ! output
586 funct_parameter(11) = zero !dpdtr! output
587 funct_parameter(12) = zero !enr! output
588
589 funct_parameter(19) = zero !dedvp
590 funct_parameter(20) = zero !pprod
591 funct_parameter(21) = zero !bthp
592 funct_parameter(22) = zero !dpdtp
593 funct_parameter(23) = zero !enp
594
595
596
597 guess = fc1 * etap/(fc1*etap+fc*etar)
598 beta = guess
599 if (beta == one) beta = fc1
600 lower_bound = epsilon(beta)
601 upper_bound = one - epsilon(beta)
602 beta = brent_algo(lower_bound,upper_bound,100*epsilon(one),jwl_eos_delta,25,funct_parameter)
603
604 etar = fc1 * etac/beta
605 CALL jwl_eos_state(ar,br,r1r,r2r,r3r,cvr,etar,tmp,dedvr,preac,bthr,dpdtr,enr)
606 etap = fc*etac/(one-beta)
607 CALL jwl_eos_state(ap,bp,r1p,r2p,r3p,cvp,etap,tmp,dedvp,pprod,bthp,dpdtp,enp)
608c pressure
609 pres = (preac+pprod)*half
610c
611c local derivatives
612 transr = bthr / beta
613 transp = bthp / (one - beta)
614 trans = transr*etar + transp*etap
615c dbdeta = d(beta)/d(eta) at equilibrium
616 dbdeta = (transr*fc1 - transp*fc) / trans
617c detar = d(etar)/d(eta) at equilibrium
618 detar = (fc1 - etar*dbdeta) / beta
619 detap = (fc + etap*dbdeta) / (one - beta)
620c dedr = d(er)/d(etar)
621 dedr = -dedvr / etar**2
622 dedp = -dedvp / etap**2
623c bth = d(p)/d(eta) at constant temperature
624 bth = (bthr*detar + bthp*detap)*half
625c dbdt = d(beta)/d(tmp) at equilibrium
626 dbdt = (dpdtr-dpdtp) / trans
627 detrdt = -etar*dbdt/beta
628 detpdt = etap*dbdt/(one-beta)
629 dpdt = (dpdtr+dpdtp+dbdt*(transp*etap-transr*etar))*half
630c dbdf = d(beta)/d(fc) at equilibrium
631 dbdf = -(transr+transp)*etac/trans
632c detrdf = d(etar)/d(fc) at equilibrium
633 detrdf = -(etac+etar*dbdf) / beta
634 detpdf = (etac+etap*dbdf) / (one-beta)
635c
636c output derivatives
637c dedv = de/dv (for constant T & fc values)
638 dedv = -etac**2 * (fc1*dedr*detar+fc*dedp*detap)
639c cv = de/dT(for constant V & fc values)
640 cvr0 = cvr + dedr*detrdt
641 cvp0 = cvp + dedp*detpdt
642 cv = fc1*cvr0 + fc*cvp0 !cv mixture
643c heat = -de/d(fc) (for constant T & V values)
644 heat = enq-enp+enr - fc*dedp*detpdf-fc1*dedr*detrdf
645c dpde = dp/de (for constant V & fc values)
646 dpde = dpdt / cv
647c dpdmu = dp/d(eta) (for constant E & fc values)
648 dpdmu = bth + dpdt*dedv/cv/etac**2
649C-----------------------------------------------
650 RETURN
651 END SUBROUTINE mixture_equilibrium
652
my_real
#define my_real
Definition cppsort.cpp:32
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
burn_fraction_ireac2
subroutine burn_fraction_ireac2(i_, b_, x_, g1, d_, y_, g2, z_, g_, c_, e_, figmax, fg1max, fg2min, cappa, chi, ccrit, tol, artvisc, dt, pres, etar, fc, cburn, oldfc)
Definition sigeps41.F:423
sigeps41
subroutine sigeps41(output, nel, nuparam, nuvar, time, timestep, uparam, rho0, rho, volume, eint, sigoxx, sigoyy, sigozz, signxx, signyy, signzz, signxy, signyz, signzx, sigvxx, sigvyy, sigvzz, sigvxy, sigvyz, sigvzx, soundsp, viscmax, uvar, deltvol, bfrac, temp)
Definition sigeps41.F:44
burn_fraction_ireac1
subroutine burn_fraction_ireac1(i_, b_, x_, g1, d_, y_, cappa, artvisc, dt, pres, etar, fc, cburn, oldfc)
Definition sigeps41.F:359
mixture_equilibrium
subroutine mixture_equilibrium(ar, br, r1r, r2r, r3r, cvr, etar, ap, bp, r1p, r2p, r3p, cvp, etap, dpdmu, dpde, ftol, enq, tmp, heat, vol_rel, fc, cv, dedv, pres, beta, funct_parameter)
Definition sigeps41.F:507