OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
jwl51.F File Reference
#include "implicit_f.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine jwlun51 (time, xl, tburn, uparam, dd, mu, mup1, volume, dvol, v1old, eint1, viscmax, q1, pext, p1, pm1, rho, rho10, mas1, ssp1, qa, qb, bfrac, dpde)
subroutine jwl51 (uparam, v1, v1old, mu1, mup1, eint1, pext, p1, pm1, rho1, rho10, mas1, ssp1, bfrac, v10, dpde1)

Function/Subroutine Documentation

◆ jwl51()

subroutine jwl51 ( dimension(*), intent(in) uparam,
intent(in) v1,
intent(in) v1old,
intent(inout) mu1,
mup1,
intent(inout) eint1,
intent(in) pext,
intent(inout) p1,
intent(in) pm1,
intent(inout) rho1,
intent(in) rho10,
intent(in) mas1,
intent(inout) ssp1,
intent(in) bfrac,
intent(in) v10,
intent(inout) dpde1 )

Definition at line 144 of file jwl51.F.

148C-----------------------------------------------
149C I m p l i c i t T y p e s
150C-----------------------------------------------
151#include "implicit_f.inc"
152C-----------------------------------------------
153C I N P U T O U T P U T A r g u m e n t s
154C-----------------------------------------------
155 my_real,INTENT(IN) :: v1,v1old,
156 . pext,pm1,
157 . rho10,mas1,
158 . uparam(*),bfrac, v10
159 my_real,INTENT(INOUT) :: rho1, mu1, eint1, p1, ssp1, dpde1
160C-----------------------------------------------
161C L o c a l V a r i a b l e s
162C-----------------------------------------------
163 INTEGER IBFRAC
164 my_real aa,bb,p0,vdet,bhe,b1,b2,w1,r1,r2,r1m,er1m,r2m,er2m,
165 . mup1,c11,c01,ssp_products,ssp_unreacted,
166 . psol, pgas, psol_min, pgas_min,dpdmu
167C-----------------------------------------------
168 vdet = uparam(42)
169 bhe = uparam(44)
170 b1 = uparam(45)
171 c01 = uparam(49)
172 c11 = uparam(50)
173 b2 = uparam(51)
174 r1 = uparam(52)
175 r2 = uparam(53)
176 w1 = uparam(54)
177 ibfrac = nint(uparam(68))
178C------------------------
179 rho1 = mas1/v1
180 mup1 = rho1/rho10
181 mu1 = mup1 - one
182
183 r1m = r1/mup1
184 r2m = r2/mup1
185 er1m = exp(-r1m)
186 er2m = exp(-r2m)
187
188 aa = w1*mup1/v10 !W1/V1 same digits this way
189 aa = aa
190 bb = half*(v1-v1old)
191 p0 = b1*(one-w1/r1m)*er1m + b2*(one-w1/r2m)*er2m
192 p1 = p0 - pext + aa * eint1
193
194 !--------------------------------!
195 ! Linear and jwl eos !
196 !--------------------------------!
197 psol = c01+c11*mu1 !linear eos relative pressure
198 psol_min = pm1 !p<0 allowed for solid phase. Default : -EP30
199 psol = max(psol,psol_min)
200
201 pgas = p1 !jwl eos relative to Pext
202 pgas_min = -pext !p>0 for detonation products
203 pgas = max(pgas,pgas_min)
204
205 p1 = bfrac*pgas + (one-bfrac)*psol
206
207 !--------------------------------!
208 ! Sound Speed !
209 !--------------------------------!
210 dpdmu = b1*er1m*( (-w1*mup1/r1) + r1m - w1) + b2*er2m*( (-w1*mup1/r2) + r2m - w1)
211 . + w1*eint1/v1 + (pgas+pext)*w1
212 dpdmu = abs(dpdmu) / mup1
213 ssp_products = sqrt(dpdmu/rho10)
214 ssp_unreacted = sqrt(c11/rho10)
215 ssp1 = (one-bfrac)*ssp_unreacted + bfrac*ssp_products
216
217 dpde1 = bfrac * w1/mup1
218
219 RETURN
my_real
#define my_real
Definition cppsort.cpp:32
max
#define max(a, b)
Definition macros.h:21

◆ jwlun51()

subroutine jwlun51 ( intent(in) time,
intent(in) xl,
intent(in) tburn,
dimension(*), intent(in) uparam,
intent(in) dd,
intent(inout) mu,
mup1,
intent(in) volume,
intent(inout) dvol,
intent(in) v1old,
intent(in) eint1,
intent(inout) viscmax,
intent(inout) q1,
intent(in) pext,
intent(inout) p1,
intent(in) pm1,
intent(inout) rho,
intent(in) rho10,
intent(in) mas1,
intent(inout) ssp1,
intent(in) qa,
intent(in) qb,
intent(inout) bfrac,
intent(inout) dpde )

Definition at line 28 of file jwl51.F.

32C-----------------------------------------------
33C I m p l i c i t T y p e s
34C-----------------------------------------------
35#include "implicit_f.inc"
36C-----------------------------------------------
37C I N P U T O U T P U T A r g u m e n t s
38C-----------------------------------------------
39 my_real, INTENT(IN) :: time,xl,tburn,dd,
40 . volume,v1old,eint1,
41 . pext,pm1,
42 . rho10,mas1,
43 . uparam(*),qa,qb
44 my_real, INTENT(INOUT) :: ssp1, p1, dvol, rho, bfrac, mu, viscmax, q1, dpde
45C-----------------------------------------------
46C L o c a l V a r i a b l e s
47C-----------------------------------------------
48 INTEGER IBFRAC
49 my_real aa,bb,p0,vdet,bhe,b1,b2,w1,r1,r2,r1m,er1m,r2m,er2m,
50 . qal,qbl,dpdmu,mup1,c01,c11,
51 . psol, pgas, psol_min, pgas_min, ssp_unreacted, ssp_reacted
52C-----------------------------------------------
53 vdet = uparam(42)
54 bhe = uparam(44)
55 b1 = uparam(45)
56 c01 = uparam(49)
57 c11 = uparam(50)
58 b2 = uparam(51)
59 r1 = uparam(52)
60 r2 = uparam(53)
61 w1 = uparam(54)
62 ibfrac= nint(uparam(68))
63C
64 IF(r1 == zero) r1=ep30
65 IF(r2 == zero) r2=ep30
66C
67 dvol = volume - v1old
68C
69 !--------------------------------!
70 ! Calculation of BFRAC in [0,1] !
71 !--------------------------------!
72 rho = mas1 / volume
73 IF(bfrac < one) THEN
74 bfrac = zero
75 IF(ibfrac/=1 .AND. time > -tburn) bfrac = vdet*(time+tburn)*two_third/xl
76 IF(ibfrac/=2) bfrac = max( bfrac , bhe * (one - rho10/rho) )
77 IF(bfrac < em03) THEN
78 bfrac = zero
79 ELSEIF(bfrac > one) THEN
80 bfrac = one
81 ENDIF
82 ENDIF
83
84 !--------------------------------!
85 ! SSP & ARTIFICIAL VISCO !
86 !--------------------------------!
87 mup1 = rho/rho10
88 mu = mup1 - one
89 r1m = r1/mup1
90 r2m = r2/mup1
91 er1m = exp(-r1m)
92 er2m = exp(-r2m)
93 aa = w1/volume
94 p0 = b1*(one-w1/r1m)*er1m + b2*(one-w1/r2m)*er2m
95 p1 = p0 + aa*eint1 !total jwl pressure for ssp
96 dpdmu = b1*er1m*( (-w1*mup1/r1) + r1m - w1) + b2*er2m*( (-w1*mup1/r2) + r2m - w1) + w1*eint1/volume +p1*w1
97 dpdmu = abs(dpdmu) / mup1 ! if DPDMU <0 => numerical error during energy integration (increase iteration number or reduce submaterial volume change ratio)
98 ssp_reacted = sqrt(dpdmu/rho10)
99 ssp_unreacted = sqrt(c11/rho10)
100 ssp1 = max(bfrac*ssp_reacted,(one-bfrac)*ssp_unreacted)
101 qal = qa*xl
102 qal = qal*qal
103 qbl = qb*xl
104 viscmax = rho*(qal*max(zero,dd) + qbl*ssp1)
105 q1 = viscmax*max(zero,dd)
106 bb = half*(volume-v1old)
107 aa = aa
108 p1 = ( p0-pext + aa*eint1 )! / (ONE+AA*BB)
109
110
111 !--------------------------------!
112 ! Linear and jwl eos !
113 !--------------------------------!
114 psol = c01+c11*mu !linear eos relative pressure
115 psol_min = pm1 !p<0 allowed for solid phase. Default : -EP30
116 psol = max(psol,psol_min)
117
118 pgas = p1 !jwl eos relative to Pext
119 pgas_min = -pext !p>0 for detonation products
120 pgas = max(pgas,pgas_min)
121
122 p1 = bfrac*pgas + (one-bfrac)*psol
123
124
125 !--------------------------------!
126 ! Update SSP with current state !
127 !--------------------------------!
128 dpdmu = b1*er1m*( (-w1*mup1/r1) + r1m - w1) + b2*er2m*( (-w1*mup1/r2) + r2m - w1) + w1*eint1/volume +(p1+pext)*w1
129 dpdmu = abs(dpdmu) / mup1 ! if DPDMU <0 => numerical error during energy integration (increase iteration number or reduce submaterial volume change ratio)
130 ssp_reacted = sqrt(dpdmu/rho10)
131 ssp_unreacted = sqrt(c11/rho10)
132 ssp1 = max(bfrac*ssp_reacted,(one-bfrac)*ssp_unreacted)
133 dpde = bfrac*w1/mup1
134C
135 RETURN