OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
jwl51.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!|| jwlun51 ../engine/source/materials/mat/mat051/jwl51.F
25!||--- called by ------------------------------------------------------
26!|| sigeps51 ../engine/source/materials/mat/mat051/sigeps51.F90
27!||====================================================================
28 SUBROUTINE jwlun51 (TIME,XL,TBURN,UPARAM,DD,MU,MUP1,
29 . VOLUME,DVOL,V1OLD,EINT1,VISCMAX,
30 . Q1,PEXT,P1,PM1,
31 . RHO,RHO10,MAS1,SSP1, QA,QB,BFRAC)
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
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! EINT1 = EINT1 - (P1OLD+PEXT+PEXT)*BB
108 aa = aa
109 p1 = ( p0-pext + aa*eint1 )! / (ONE+AA*BB)
110
111
112 !--------------------------------!
113 ! Linear and jwl eos !
114 !--------------------------------!
115 psol = c01+c11*mu !linear eos relative pressure
116 psol_min = pm1 !p<0 allowed for solid phase. Default : -EP30
117 psol = max(psol,psol_min)
118
119 pgas = p1 !jwl eos relative to Pext
120 pgas_min = -pext !p>0 for detonation products
121 pgas = max(pgas,pgas_min)
122
123 p1 = bfrac*pgas + (one-bfrac)*psol
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)
133C
134 RETURN
135 END
136
137
138!||====================================================================
139!|| jwl51 ../engine/source/materials/mat/mat051/jwl51.F
140!||--- called by ------------------------------------------------------
141!|| sigeps51 ../engine/source/materials/mat/mat051/sigeps51.F90
142!||====================================================================
143 SUBROUTINE jwl51 (UPARAM,
144 . V1,V1OLD,MU1,MUP1,EINT1,
145 . P1OLD,PEXT,P1,PM1,
146 . RHO1,RHO10,MAS1,SSP1,DVDP1,DPDV1,BFRAC,V10, FLAG)
147C-----------------------------------------------
148C I m p l i c i t T y p e s
149C-----------------------------------------------
150#include "implicit_f.inc"
151C-----------------------------------------------
152C I N P U T O U T P U T A r g u m e n t s
153C-----------------------------------------------
154 my_real,INTENT(IN) :: v1,v1old,
155 . p1old,pext,pm1,
156 . rho10,mas1,
157 . uparam(*),bfrac, v10
158 INTEGER,INTENT(IN) :: FLAG
159 my_real,INTENT(INOUT) :: RHO1, MU1, EINT1, P1, SSP1, DVDP1,DPDV1
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,dvdp1i,c11,c01,ssp_products,ssp_unreacted,
166 . psol, pgas, psol_min, pgas_min,dpdmu,dpdv1_reacted,dpdv1_unreacted
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 dvdp1i = dvdp1
180 rho1 = mas1/v1
181 mup1 = rho1/rho10
182 mu1 = mup1 - one
183
184 r1m = r1/mup1
185 r2m = r2/mup1
186 er1m = exp(-r1m)
187 er2m = exp(-r2m)
188
189
190 aa = w1*mup1/v10 !W1/V1 same digits this way
191 aa = aa
192 bb = half*(v1-v1old)
193 IF (flag == 1) eint1 = eint1 - (p1old+pext+pext)*bb
194 p0 = b1*(one-w1/r1m)*er1m + b2*(one-w1/r2m)*er2m
195 IF (flag == 1) THEN
196 p1 = ( p0-pext + aa*eint1 ) / (one+aa*bb)
197 ELSE
198 p1 = p0 - pext + aa * eint1
199 ENDIF
200
201 !--------------------------------!
202 ! Linear and jwl eos !
203 !--------------------------------!
204 psol = c01+c11*mu1 !linear eos relative pressure
205 psol_min = pm1 !p<0 allowed for solid phase. Default : -EP30
206 psol = max(psol,psol_min)
207
208 pgas = p1 !jwl eos relative to Pext
209 pgas_min = -pext !p>0 for detonation products
210 pgas = max(pgas,pgas_min)
211
212 p1 = bfrac*pgas + (one-bfrac)*psol
213 IF (flag == 1) eint1 = eint1 - p1*bb
214 IF (flag == 1) eint1 = max(eint1, zero)
215
216 !--------------------------------!
217 ! Sound Speed !
218 !--------------------------------!
219 dpdmu = b1*er1m*( (-w1*mup1/r1) + r1m - w1) + b2*er2m*( (-w1*mup1/r2) + r2m - w1)
220 . + w1*eint1/v1 + (pgas+pext)*w1
221 dpdmu = abs(dpdmu) / mup1
222 ssp_products = sqrt(dpdmu/rho10)
223 ssp_unreacted = sqrt(c11/rho10)
224 ssp1 = (one-bfrac)*ssp_unreacted + bfrac*ssp_products
225
226 !--------------------------------!
227 ! DPDV !
228 !--------------------------------!
229 dpdv1_reacted = -dpdmu*mup1/v1
230 dpdv1_unreacted = -c11*mup1/v1
231 dpdv1 = bfrac*dpdv1_reacted + (one-bfrac)*dpdv1_unreacted
232
233 IF(abs(dpdv1)<em20)THEN
234 dvdp1 = zero
235 ELSE
236 dvdp1 = one/dpdv1
237 ENDIF
238
239 RETURN
240 END
241
my_real
#define my_real
Definition cppsort.cpp:32
jwl51
subroutine jwl51(uparam, v1, v1old, mu1, mup1, eint1, p1old, pext, p1, pm1, rho1, rho10, mas1, ssp1, dvdp1, dpdv1, bfrac, v10, flag)
Definition jwl51.F:147
jwlun51
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)
Definition jwl51.F:32
max
#define max(a, b)
Definition macros.h:21