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,DPDE)
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
136 END
137
138
139!||====================================================================
140!|| jwl51 ../engine/source/materials/mat/mat051/jwl51.F
141!||--- called by ------------------------------------------------------
142!|| sigeps51 ../engine/source/materials/mat/mat051/sigeps51.F90
143!||====================================================================
144 SUBROUTINE jwl51 (UPARAM,
145 . V1,V1OLD,MU1,MUP1,EINT1,
146 . PEXT,P1,PM1,
147 . RHO1,RHO10,MAS1,SSP1,BFRAC,V10, DPDE1)
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
220 END
221
my_real
#define my_real
Definition cppsort.cpp:32
jwl51
subroutine jwl51(uparam, v1, v1old, mu1, mup1, eint1, pext, p1, pm1, rho1, rho10, mas1, ssp1, bfrac, v10, dpde1)
Definition jwl51.F:148
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, dpde)
Definition jwl51.F:32
max
#define max(a, b)
Definition macros.h:21