OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
multi_submatlaw.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!|| multi_submatlaw_mod ../engine/source/multifluid/multi_submatlaw.F
25!||--- called by ------------------------------------------------------
26!|| multi_inlet_ebcs ../engine/source/multifluid/multi_inlet_ebcs.F
27!|| multi_muscl_compute_pressure ../engine/source/multifluid/multi_muscl_compute_pressure.F90
28!|| multi_muscl_fluxes_computation ../engine/source/multifluid/multi_muscl_fluxes_computation.F
29!|| multi_nrf_ebcs ../engine/source/multifluid/multi_nrf_ebcs.F
30!|| multi_pressure_equilibrium ../engine/source/multifluid/multi_pressure_equilibrium.F
31!|| multi_propellant_ebcs ../engine/source/multifluid/multi_propellant_ebcs.F90
32!|| multi_solve_eint ../engine/source/multifluid/multi_solve_eint.F90
33!||====================================================================
34 MODULE multi_submatlaw_mod
35 CONTAINS
36!||====================================================================
37!|| multi_submatlaw ../engine/source/multifluid/multi_submatlaw.F
38!||--- called by ------------------------------------------------------
39!|| multi_inlet_ebcs ../engine/source/multifluid/multi_inlet_ebcs.F
40!|| multi_muscl_compute_pressure ../engine/source/multifluid/multi_muscl_compute_pressure.F90
41!|| multi_nrf_ebcs ../engine/source/multifluid/multi_nrf_ebcs.F
42!|| multi_pressure_equilibrium ../engine/source/multifluid/multi_pressure_equilibrium.f
43!|| multi_propellant_ebcs ../engine/source/multifluid/multi_propellant_ebcs.F90
44!|| multi_solve_eint ../engine/source/multifluid/multi_solve_eint.F90
45!||--- calls -----------------------------------------------------
46!|| arret ../engine/source/system/arret.F
47!|| eosmain ../common_source/eos/eosmain.F
48!||--- uses -----------------------------------------------------
49!|| alemuscl_mod ../common_source/modules/ale/alemuscl_mod.F
50!|| eosmain_mod ../common_source/eos/eosmain.F
51!|| matparam_def_mod ../common_source/modules/mat_elem/matparam_def_mod.F90
52!||====================================================================
53 SUBROUTINE multi_submatlaw(
54 1 IFLAG, MATLAW, LOCAL_MATID, NEL,
55 2 EINT, PRES, RHO, SSP,
56 3 VOL, GRUN, PM, IPM,
57 4 NPROPM, NPROPMI, BUFMAT, OFF,
58 5 THETA, BURNFRAC, BURNTIME, DELTAX,
59 6 CURRENT_TIME,SIGOLD, SNPC, STF,
60 7 NPF, TF, VAREOS, NVAREOS,
61 8 MAT_PARAM, NVARTMP_EOS, VARTMP_EOS, NUMMAT,
62 9 ABURN)
63C-----------------------------------------------
64C M o d u l e s
65C-----------------------------------------------
66 USE matparam_def_mod , ONLY : matparam_struct_
67 USE eosmain_mod , ONLY : eosmain
68 USE alemuscl_mod , only : alemuscl_param
69C-----------------------------------------------
70C I m p l i c i t T y p e s
71C-----------------------------------------------
72#include "implicit_f.inc"
73C-----------------------------------------------
74C C o m m o n B l o c k s
75C-----------------------------------------------
76#include "com08_c.inc"
77C-----------------------------------------------
78C D u m m y A r g u m e n t s
79C-----------------------------------------------
80 INTEGER,INTENT(IN) :: SNPC, STF,NUMMAT
81 INTEGER, INTENT(IN) :: IFLAG
82 INTEGER, INTENT(IN) :: MATLAW, LOCAL_MATID, NEL, NPROPM, NPROPMI
83 my_real, INTENT(IN) :: PM(NPROPM, NUMMAT),SIGOLD(NEL,6)
84 INTEGER, INTENT(IN) :: IPM(NPROPMI, NUMMAT)
85 my_real, INTENT(IN) :: VOL(NEL), OFF(NEL)
86 my_real, INTENT(OUT) :: GRUN(NEL)
87 my_real, INTENT(INOUT) :: EINT(NEL), RHO(NEL), THETA(NEL)
88 my_real, INTENT(OUT) :: PRES(NEL), SSP(NEL)
89 my_real, INTENT(IN) :: burntime(nel), deltax(nel), current_time
90 my_real, INTENT(INOUT) :: burnfrac(nel), bufmat(*)
91 INTEGER, INTENT(IN) :: NPF(SNPC),NVAREOS
92 my_real, INTENT(IN) :: tf(stf),vareos(nvareos*nel)
93 TYPE(matparam_struct_), INTENT(IN) :: MAT_PARAM !material data structure
94 INTEGER,INTENT(IN) :: NVARTMP_EOS
95 INTEGER,INTENT(INOUT) :: VARTMP_EOS(NEL,NVARTMP_EOS)
96 my_real,INTENT(INOUT) :: aburn(nel) !< after burning
97C-----------------------------------------------
98C L o c a l V a r i a b l e s
99C-----------------------------------------------
100 INTEGER :: II, EOSTYPE
101 my_real :: c0, c1, rho0(nel)
102 my_real :: mu(nel), mu2(nel), espe(nel), dpde(nel), ecold(nel), dvol(nel),
103 . df(nel), psh(nel),muold(nel)
104 INTEGER :: MAT(NEL)
105 my_real :: gamma, pstar
106 my_real :: b1, b2, r1, r2, w1, vdet, bhe
107 my_real :: r1v, r2v, wdr1v, wdr2v, dr1v, er1v, er2v
108 my_real :: tb, bfrac1, bfrac2
109 INTEGER :: IBFRAC
110 INTEGER :: LFT, LLT
111 my_real :: qopt,eadd,tbegin,tend,rr,a,m,n,rr2,alpha_unit,lambda !< jwl afterburning
112 my_real :: pold(nel)
113 my_real :: fac_pred
114C-----------------------------------------------
115C B e g i n n i n g o f s u b r o u t i n e
116C-----------------------------------------------
117 fac_pred = one
118 IF(alemuscl_param%IALEMUSCL /= 0) fac_pred = half ! with prediction /correction steps, Miller's extension for JWL must use a half time step
119 lft = 1
120 llt = nel
121 rho0(1:nel) = pm(1, local_matid)
122 mat(1:nel) = local_matid
123 dvol(1:nel) = zero
124 mu2(1:nel) = zero
125 mu(1:nel) = zero
126 espe(1:nel) = zero
127 dpde(1:nel) = zero
128 ecold(1:nel) = zero
129 dvol(1:nel) = zero
130 df(1:nel) = one
131 psh(1:nel) = zero
132 muold(1:nel) = zero
133 gamma = zero
134 pstar = zero
135 DO ii = 1, nel
136 IF (vol(ii) > zero) THEN
137 mu(ii) = rho(ii) / rho0(ii) - one
138 df(ii) = rho0(ii) / rho(ii)
139 psh(ii) = pm(88, local_matid)
140 ELSE
141 mu(ii) = zero
142 df(ii) = one
143 psh(ii) = -1e20
144 ENDIF
145 ENDDO
146 SELECT CASE (matlaw)
147 CASE (5)
148C JWL material
149 b1 = pm(33, local_matid)
150 b2 = pm(34, local_matid)
151 r1 = pm(35, local_matid)
152 r2 = pm(36, local_matid)
153 w1 = pm(45, local_matid)
154 vdet = pm(38, local_matid)
155 bhe = pm(40, local_matid)
156 c0 = pm(43, local_matid)
157 c1 = pm(44, local_matid)
158 qopt = pm(42, local_matid)
159 eadd = pm(160,local_matid)
160 tbegin = pm(161,local_matid)
161 tend = pm(162,local_matid)
162 rr = pm(163,local_matid)
163 a = pm(164,local_matid)
164 m = pm(165,local_matid)
165 n = pm(166,local_matid)
166 rr2 = pm(167,local_matid)
167 alpha_unit = pm(168,mat(1))
168 ibfrac = nint(pm(41, local_matid))
169 grun(1:nel) = zero
170 DO ii = 1, nel
171 IF (vol(ii) > zero) THEN
172 !================ BURNT FRACTION ============================!
173 IF (iflag == 1) THEN
174 IF (burnfrac(ii) < one) THEN
175 tb = - burntime(ii)
176 bfrac1 = zero
177 bfrac2 = zero
178 IF (ibfrac /= 1 .AND. current_time > tb) THEN
179 bfrac1 = one
180 IF (deltax(ii) > zero) THEN
181 bfrac1 = vdet * (current_time - tb) / three_half / deltax(ii)
182 ENDIF
183 ENDIF
184 IF (ibfrac /= 2) THEN
185 bfrac2 = bhe * (one - df(ii))
186 ENDIF
187 !BURNFRAC(II) = MAX(BFRAC1, BFRAC2, BURNFRAC(II))
188 burnfrac(ii) = max(bfrac1, bfrac2)
189 IF (burnfrac(ii) < em04) THEN
190 burnfrac(ii) = zero
191 ENDIF
192 IF (burnfrac(ii) > one) THEN
193 burnfrac(ii) = one
194 ENDIF
195 ENDIF
196 ENDIF
197 !================ END BURNT FRACTION ============================!
198 ENDIF
199 ENDDO
200
201 ! POLD
202 pold(:) = zero
203 DO ii = 1, nel
204 IF (vol(ii) > zero) THEN
205 r1v = b1 * w1 / (r1 * df(ii))
206 r2v = b2 * w1 / (r2 * df(ii))
207 wdr1v = b1 - r1v
208 wdr2v = b2 - r2v
209 dr1v = w1 * eint(ii)
210 er1v = exp(-r1 * df(ii))
211 er2v = exp(-r2 * df(ii))
212 pold(ii) = - psh(ii) + wdr1v * er1v + wdr2v * er2v + dr1v
213 pold(ii) = burnfrac(ii) * pold(ii) + (one - burnfrac(ii)) * (-psh(ii)+(c0 + c1 * mu(ii)))
214 ENDIF
215 ENDDO
216
217 !================ AFTERBURNING ============================!
218 IF(eadd /= zero)THEN
219 SELECT CASE (nint(qopt))
220 CASE(0) !=== instantaneous release
221 DO ii=1,nel
222 IF (vol(ii) > zero) THEN
223 lambda = zero
224 IF(current_time > tend .AND. aburn(ii)==one)THEN
225 lambda = one
226 aburn(ii) = one
227 ELSEIF (current_time <= tbegin)THEN
228 lambda = zero
229 aburn(ii) = zero
230 ELSE
231 lambda = one
232 eint(ii) = eint(ii)+(lambda-aburn(ii))*eadd*max(em20,vol(ii)/df(ii))
233 aburn(ii) = one
234 ENDIF
235 ENDIF !VOL(II) > ZERO
236 ENDDO
237 CASE(1) !=== afterburning with constant rate from Tbegin to Tend
238 DO ii=1,nel
239 IF (vol(ii) > zero) THEN
240 lambda = zero
241 IF(current_time > tend .AND. aburn(ii)==one)THEN
242 lambda = one
243 aburn(ii) = one
244 ELSEIF (current_time <= tbegin)THEN
245 lambda = zero
246 aburn(ii) = zero
247 ELSE
248 lambda = (current_time-tbegin)*rr
249 lambda = min(one,lambda)
250 eint(ii) = eint(ii)+(lambda-aburn(ii))*eadd*max(em20,vol(ii)/df(ii))
251 aburn(ii) = lambda
252 ENDIF
253 ENDIF !VOL(II) > ZERO
254 ENDDO
255 CASE(2) !=== afterburning with linear rate from Tbegin to Tend
256 DO ii=1,nel
257 IF (vol(ii) > zero) THEN
258 lambda = zero
259 IF(current_time > tend .AND. aburn(ii)==one)THEN ! .AND. ABURN(II)==ONE needed to add last increment
260 lambda = one
261 aburn(ii) = one
262 ELSEIF (current_time <= tbegin)THEN
263 lambda = zero
264 aburn(ii) = zero
265 ELSE
266 lambda = half*rr*current_time**2 - rr*tbegin*current_time + rr2
267 lambda = max(zero,min(one,lambda))
268 eint(ii) = eint(ii)+(lambda-aburn(ii))*eadd*max(em20,vol(ii)/df(ii))
269 aburn(ii) = lambda
270 ENDIF
271 ENDIF !VOL(II) > ZERO
272 ENDDO
273 CASE(3) !=== Miller s extension, rate is depedent on Pressure
274 DO ii=1,nel
275 IF (vol(ii) > zero) THEN
276 lambda = zero
277 IF(pold(ii)+psh(ii) > zero )THEN
278 lambda=aburn(ii)+ fac_pred*dt1*a*exp( m*log(one+aburn(ii)) )*exp(n*log(alpha_unit*(pold(ii)+psh(ii))))
279 lambda = max(lambda,zero)
280 lambda = min(lambda,one)
281 eint(ii) = eint(ii)+(lambda-aburn(ii))*eadd*max(em20,vol(ii)/df(ii))/fac_pred
282 aburn(ii)= lambda
283 ENDIF
284 ENDIF !VOL(II) > ZERO
285 ENDDO
286 END SELECT
287 ENDIF !EADD /= ZERO
288 !================ END AFTERBURNING ============================!
289
290 DO ii = 1, nel
291 IF (vol(ii) > zero) THEN
292 r1v = b1 * w1 / (r1 * df(ii))
293 r2v = b2 * w1 / (r2 * df(ii))
294 wdr1v = b1 - r1v
295 wdr2v = b2 - r2v
296 dr1v = w1 * eint(ii)
297 er1v = exp(-r1 * df(ii))
298 er2v = exp(-r2 * df(ii))
299 pres(ii) = - psh(ii) + wdr1v * er1v + wdr2v * er2v + dr1v
300 ssp(ii) = b1 * er1v * (-w1 / df(ii) / r1 + r1 * df(ii) - w1)
301 . + b2 * er2v * (-w1 / df(ii) / r2 + r2 * df(ii) - w1)
302 . + dr1v + (pres(ii) + psh(ii)) * w1
303 pres(ii) = burnfrac(ii) * pres(ii) +
304 . (one - burnfrac(ii)) * (-psh(ii)+(c0 + c1 * mu(ii)))
305 grun(ii) = burnfrac(ii) * w1
306 ssp(ii) = burnfrac(ii) * (ssp(ii) * df(ii) / rho0(ii)) +
307 . (one - burnfrac(ii)) * (c1 / rho0(ii))
308 !SSP(II) = SQRT(SSP(II))
309 !SSP(II) = MAX(SSP(II), VDET * (ONE - BURNFRAC(II)))
310 ENDIF
311 ENDDO
312 CASE (3, 4, 6, 49)
313 eostype = ipm(4, local_matid)
314 muold(1:nel) = mu(1:nel)
315C -> FLAG = 2 : compute pressure as a function of rho and e, as well as dpdm and dpde
316 CALL eosmain(2 , nel , eostype , pm , off , eint,
317 . rho , rho0 , mu , mu2 , espe ,
318 . dvol , df , vol , mat , psh ,
319 . pres , ssp , dpde , theta ,
320 . bufmat , sigold , muold , matlaw ,
321 . npf , tf , vareos , nvareos, mat_param,
322 . burnfrac, nvartmp_eos, vartmp_eos)
323C SQUARE sound speed : c^2 = 1/rho0 * dp/dmu
324 DO ii = 1, nel
325 IF (vol(ii) > zero) THEN
326 ssp(ii) = ssp(ii) / rho0(ii)
327 grun(ii) = dpde(ii) / (one + mu(ii))
328 ELSE
329 ssp(ii) = zero
330 grun(ii) = zero
331 ENDIF
332 ENDDO
333 CASE DEFAULT
334 print*, "LAW", matlaw, "NOT COMPATIBLE WITH LAW 151 (MULTIFLUID)"
335 CALL arret(2)
336 END SELECT
337C-----------------------------------------------
338C E n d o f s u b r o u t i n e
339C-----------------------------------------------
340 END SUBROUTINE multi_submatlaw
341
342
343 END MODULE multi_submatlaw_mod
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
multi_pressure_equilibrium
subroutine multi_pressure_equilibrium(timestep, elbuf_tab, iparg, itask, ixs, ixq, ixtg, pm, ipm, multi_fvm, current_time, bufmat, npf, tf, nummat, matparam)
Definition multi_pressure_equilibrium.F:38
alemuscl_mod
Definition alemuscl_mod.F:54
alemuscl_mod::alemuscl_param
type(alemuscl_param_) alemuscl_param
Definition alemuscl_mod.F:75
eosmain_mod
Definition eosmain.F:33
eosmain_mod::eosmain
subroutine eosmain(iflag, nel, eostyp, pm, off, eint, rho, rho0, mu, mu2, espe, dvol, df, vnew, mat, psh, pnew, dpdm, dpde, theta, bufmat, sig, mu_bak, mlw, npf, tf, vareos, nvareos, mat_param, bfrac, nvartmp, vartmp)
Definition eosmain.F:80
multi_submatlaw_mod
Definition multi_submatlaw.F:34
multi_submatlaw_mod::multi_submatlaw
subroutine multi_submatlaw(iflag, matlaw, local_matid, nel, eint, pres, rho, ssp, vol, grun, pm, ipm, npropm, npropmi, bufmat, off, theta, burnfrac, burntime, deltax, current_time, sigold, snpc, stf, npf, tf, vareos, nvareos, mat_param, nvartmp_eos, vartmp_eos, nummat, aburn)
Definition multi_submatlaw.F:63
arret
subroutine arret(nn)
Definition arret.F:87