OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
m51init.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!|| m51init ../starter/source/materials/mat/mat051/m51init.f
25!||--- called by ------------------------------------------------------
26!|| matini ../starter/source/materials/mat_share/matini.F
27!||--- calls -----------------------------------------------------
28!|| m5in2 ../starter/source/initial_conditions/detonation/m5in2.F
29!|| m5in3 ../starter/source/initial_conditions/detonation/m5in3.F
30!|| nrf51ini ../starter/source/materials/mat/mat051/nrf51ini.F
31!||--- uses -----------------------------------------------------
32!|| detonators_mod ../starter/share/modules1/detonators_mod.F
33!||====================================================================
34 SUBROUTINE m51init(
35 . IPM ,DETONATORS ,PM ,TB ,
36 . NUVAR ,UVAR ,UPARAM ,X ,
37 . MAT ,IPARG ,IFORM ,IX ,NIX ,
38 . ALE_CONNECTIVITY ,BUFMAT ,RHO0 ,
39 . GBUF ,NEL ,SIG ,MAT_PARAM , NPF, TF )
40C-----------------------------------------------
41C M o d u l e s
42C-----------------------------------------------
43 USE elbufdef_mod
44 USE detonators_mod
45 USE ale_connectivity_mod
46 USE multimat_param_mod , ONLY : m51_n0phas, m51_nvphas
47 USE matparam_def_mod , ONLY : matparam_struct_
48 USE eosmain_mod , ONLY : eosmain
49C-----------------------------------------------
50C I m p l i c i t T y p e s
51C-----------------------------------------------
52#include "implicit_f.inc"
53C-----------------------------------------------
54C C o m m o n B l o c k s
55C-----------------------------------------------
56#include "com01_c.inc"
57#include "com04_c.inc"
58#include "param_c.inc"
59C-----------------------------------------------
60C D u m m y A r g u m e n t s
61C-----------------------------------------------
62 INTEGER :: IPM(NPROPMI,NUMMAT),MAT(NEL), IPARG(NPARG),IFORM,NIX,IX(NIX,*)
63 my_real :: PM(NPROPM,NUMMAT),UPARAM(*), X(3,NUMNOD), BUFMAT(*), RHO0
64 my_real , TARGET :: UVAR(NEL,NUVAR)
65 INTEGER,INTENT(IN) :: NEL
66 my_real,INTENT(INOUT) :: SIG(NEL,6)
67 TYPE(g_bufel_), INTENT(INOUT),TARGET :: GBUF
68 TYPE(detonators_struct_) :: DETONATORS
69 TYPE(t_ale_connectivity), INTENT(INOUT) :: ALE_CONNECTIVITY
70 TYPE (MATPARAM_STRUCT_), DIMENSION(NUMMAT), INTENT(INOUT) :: MAT_PARAM
71 INTEGER,INTENT(IN) :: NPF(*)
72 my_real,INTENT(IN) :: tf(*)
73C-----------------------------------------------
74C L o c a l V a r i a b l e s
75C-----------------------------------------------
76 INTEGER :: I,NUVAR, ISFLUID
77 INTEGER :: NPH,IFLG,NV46
78 INTEGER :: IMAT,MID
79 INTEGER :: KK,NUMEL
80 INTEGER :: IDX_YIELD(4)
81 INTEGER :: EOSTYPE, MATLAW
82 INTEGER :: NBMAT
83 my_real :: gg1, gg2, gg3, gg4
84 my_real :: vold,p0_(4),t0_(4),pmin_(4),e0_(4),rho0_(4),vf(4),pext,pres_0,ssp_(4),gg_(4),dpdm_(4)
85 my_real :: c0_(4),c1_(4),c2_(4),c3_(4),c4_(4),c5_(4)
86 my_real,intent(inout) :: tb(nel)
87 my_real,DIMENSION(1) :: off, eint, rho, rho_0, mu, mu2, espe, muold, burnfrac
88 my_real,DIMENSION(1) :: dvol, df, vol, pres, ssp, dpde, dpdm, theta, psh
89 INTEGER :: MAT_(1), ISUB
90 my_real :: sph1, sph2, sph3, sph4
91 my_real :: v1, v2, v3, v4, k_mix, ssp_mix, rho0_mix, e0_mix
92 my_real,DIMENSION(NEL,6) :: sigold
93 my_real :: pmin, denom
94 LOGICAL :: IS_IFORM12
95 integer, parameter :: nvartmp = 3
96 integer, parameter :: nvareos = 6
97 my_real,DIMENSION(1,nvareos) :: vareos ! pointer to UVAR
98 integer :: vartmp(1, nvartmp)
99C-----------------------------------------------
100
101 iflg = nint(uparam(55))
102 iform = nint(uparam(31))
103 sph1 = uparam(112)
104 sph2 = uparam(162)
105 sph3 = uparam(212)
106 sph4 = uparam(262)
107
108 idx_yield(1:4) = (/100,150,200,250/)
109
110 !need to initialize element buffer with initial density for this specific case.
111 is_iform12 = .false.
112 IF( nint(uparam(76)) == 12) is_iform12 = .true.
113 IF(iform == 12)THEN
114 uparam(31) = 1
115 iform = 1
116 ENDIF
117
118 !---------------------------------!
119 ! Test if only fluids are defined !
120 !---------------------------------!
121 gg1 = uparam(101)
122 gg2 = uparam(151)
123 gg3 = uparam(201)
124 gg4 = zero
125 gg_(1:4) = (/gg1,gg2,gg3,gg4/)
126 isfluid = 0
127 IF (gg1 == zero .AND. gg2 == zero .AND. gg3 == zero) isfluid=1
128 IF (isfluid==1) THEN
129 iparg(15) = 1 ! FLAG REZONE SIG => KEPT TO 1 FOR BACKWARD COMPATIBILITY
130 iparg(16) = 1 ! FLAG REZONE EPS PLAST. => KEPT TO 1 FOR BACKWARD COMPATIBILITY
131 iparg(63) = 1 ! FLAG FOR FLUID MATERIAL
132 iparg(64) = 0 ! FLAG FOR SILENT BOUNDARY FORMULATION
133 ENDIF
134
135 !---------------------------------!
136 ! Test if NRF formulation !
137 !---------------------------------!
138 IF(iform == 2 .OR. iform == 3 .OR. iform == 4 .OR. iform==5 .OR. iform == 6) THEN
139 iparg(15) = 0 ! FLAG REZONE SIG
140 iparg(16) = 0 ! FLAG REZONE EPS PLAST.
141 iparg(63) = 1 ! FLAG FOR FLUID MATERIAL
142 iparg(64) = 1 ! FLAG FOR SILENT BOUNDARY FORMULATION
143 ENDIF
144
145 !---------------------------------!
146 ! compute burning time !
147 !---------------------------------!
148 IF(iflg == 1 .AND. iparg(64) == 0)THEN !unplug for NRF
149 nph = 1
150 IF(n2d == 0)THEN
151 CALL m5in3 (pm,mat,0,detonators,tb,iparg,x,ix,nix)
152 ELSE
153 CALL m5in2 (pm,mat,0,detonators,tb,x,ix,nix)
154 ENDIF
155 ENDIF
156
157 !C==========================================
158 !C USER VARIABLES INITIALIZATION T=0
159 !C==========================================
160 c0_(1:4) = (/uparam(35:37),uparam(49)/)
161 c1_(1:4) = (/uparam(12:14),uparam(50)/)
162 c2_(1:4) = (/uparam(15:17),zero/)
163 c3_(1:4) = (/uparam(18),uparam(20:21),zero/)
164 c4_(1:4) = (/uparam(22:24),zero/)
165 c5_(1:4) = (/uparam(25:27),zero/)
166 e0_(1:4) = (/uparam(32:34),uparam(48)/)
167 t0_(1:4) = (/uparam(113),uparam(163),uparam(213),uparam(263)/)
168 rho0_(1:4) = (/uparam(09:11),uparam(47)/)
169 vf(1:4) = (/uparam(04:06),uparam(46)/)
170 pext = uparam(8)
171 p0_(1:4) = (/uparam(57:60)/)
172 pmin_(1:4) = (/uparam(39),uparam(40),uparam(41),uparam(56)/)
173 t0_(1:4) = (/uparam(idx_yield(1)+13),uparam(idx_yield(1)+13),uparam(idx_yield(1)+13),uparam(idx_yield(1)+13)/)
174
175 ! INITIAL SOUND SPEED (mu=0,E=E0)
176
177 IF(.NOT. is_iform12)THEN
178
179 ! old way with embedded polynomial EoS (backward compatibility)
180 p0_(1) = c0_(1)+c4_(1)*e0_(1)
181 p0_(2) = c0_(2)+c4_(2)*e0_(2)
182 p0_(3) = c0_(3)+c4_(3)*e0_(3)
183 p0_(4) = c0_(4)
184 dpdm_(1) = c1_(1)+c5_(1)*e0_(1)+p0_(1)*c4_(1)
185 dpdm_(2) = c1_(2)+c5_(2)*e0_(2)+p0_(2)*c4_(2)
186 dpdm_(3) = c1_(3)+c5_(3)*e0_(3)+p0_(3)*c4_(3)
187 dpdm_(4) = c1_(4)+c5_(4)*e0_(4)+p0_(4)*c4_(4)
188 ssp_(1:4)=zero
189 IF(rho0_(1) > zero)ssp_(1) = sqrt( (abs(dpdm_(1))+two*two_third*gg1)/rho0_(1) )
190 IF(rho0_(2) > zero)ssp_(2) = sqrt( (abs(dpdm_(2))+two*two_third*gg2)/rho0_(2) )
191 IF(rho0_(3) > zero)ssp_(3) = sqrt( (abs(dpdm_(3))+two*two_third*gg3)/rho0_(3) )
192 IF(rho0_(4) > zero)ssp_(4) = sqrt( (abs(dpdm_(4))+two*two_third*gg4)/rho0_(4) )
193
194 ELSE
195
196 nbmat = mat_param(mat(1))%MULTIMAT%NB
197 vareos(1,1:6) = zero
198 ssp_(1:4)=zero
199 DO imat=1,nbmat
200 mid = mat_param(mat(1))%MULTIMAT%MID(imat)
201 isub = nint(uparam(276+imat)) !convert IMAT in 1..3 -> 1..4
202 IF(mid == 0) cycle
203 off(1) = zero
204 mat_(1) = mid
205 eint(1) = e0_(isub)
206 rho(1) = rho0_(isub)
207 rho_0(1) = rho0_(isub)
208 mu(1) = zero
209 mu2(1) = zero
210 espe(1) = e0_(isub)
211 dvol(1) = zero
212 df(1) = one
213 vol(1) = one !iflag=2 in EOSMAIN : vnew > 0.
214 pres(1) = zero
215 ssp(1) = zero
216 dpde(1) = zero
217 dpdm(1) = zero
218 theta(1)= t0_(isub)
219 sigold(1,1:3) = - mat_param(mat(1))%MULTIMAT%pEOS(imat)%EOS%P0 !ideal gas vt only
220 sigold(1,4:6) = zero
221 muold(1) = zero
222 burnfrac(1) = zero
223 psh(1) = mat_param(mat(1))%MULTIMAT%pEOS(imat)%EOS%PSH
224 matlaw = ipm(2,mid)
225 pmin = pmin_(isub)
226 eostype = mat_param(mat(1))%MULTIMAT%pEOS(imat)%EOS%EOSTYPE
227 if(matlaw == 5) eostype = 15
228 vartmp(1,1:nvartmp) = 1
229 CALL eosmain(2 , 1 , eostype , pm , off , eint,
230 . rho , rho_0 , mu , mu2 , espe ,
231 . dvol , df , vol , mat_ , psh ,
232 . pres , dpdm , dpde , theta ,
233 . bufmat , sigold , muold , matlaw ,
234 . npf , tf , vareos , 6, mat_param(mid),
235 . burnfrac, nvartmp , vartmp)
236 IF(matlaw ==5) dpdm(1) = max(pm(44,mid) , dpdm(1))
237
238 IF(rho0_(isub) > zero)ssp_(isub) = sqrt( abs(dpdm(1))+two*two_third*gg_(isub)/rho0_(isub) )
239 kk = m51_n0phas + (isub-1)*m51_nvphas
240 DO i=1,nel
241 uvar(i,kk+24:kk-1+nvareos) = vareos(1, 1:6)
242 ENDDO
243 ENDDO ! next IMAT
244
245 ENDIF
246
247
248 !===========================!
249 ! material IMAT : TIME=0 !
250 ! default : !
251 ! if no /INIGRAV option !
252 !===========================!
253 DO imat=1,4
254 DO i=1,nel
255 vold = gbuf%VOL(i)*vf(imat)
256 kk = m51_n0phas + (imat-1)*m51_nvphas
257 uvar(i,1+kk) = vf(imat)
258 uvar(i,2+kk) = zero
259 uvar(i,3+kk) = zero
260 uvar(i,4+kk) = zero
261 uvar(i,5+kk) = zero
262 uvar(i,6+kk) = zero
263 uvar(i,7+kk) = zero
264 uvar(i,8+kk) = e0_(imat) ! rho.e for ACONVE
265 uvar(i,9+kk) = rho0_(imat) ! rho for ACONVE
266 uvar(i,10+kk) = zero
267 uvar(i,11+kk) = vold
268 uvar(i,12+kk) = rho0_(imat) ! rho_old IN rho OUT
269 uvar(i,14+kk) = ssp_(imat)
270 uvar(i,15+kk) = zero ! Plastic EPS
271 uvar(i,16+kk) = t0_(imat) ! temperature
272 uvar(i,17+kk) = zero ! diffuse heat
273 uvar(i,18+kk) = p0_(imat)
274 uvar(i,19+kk) = zero
275 uvar(i,20+kk) = rho0_(imat) !rho0 initial state is element dependent and can be updated by /INIGRAV
276 uvar(i,21+kk) = e0_(imat) !E0
277 uvar(i,22+kk) = zero !SSP0 to introduce if needed
278 uvar(i,23+kk) = vf(imat) !alpha0
279 !24 : vareos : defined above
280 ENDDO
281 ENDDO
282 !particular case of submaterial 4 (burning time to initialize)
283 imat = 4
284 kk = m51_n0phas+(imat-1)*m51_nvphas
285 uvar(1:nel,15+kk) = uvar(1:nel,1) ! Tdet
286 !
287 DO i=1,nel
288 pres_0 = vf(1)*p0_(1)+vf(2)*p0_(2)+vf(3)*p0_(3)+vf(4)*p0_(4)
289 uvar(i,4) = pres_0
290 pm(31,mat(i)) = pres_0
291 pm(104,mat(i)) = pres_0
292 sig(i,1) = -pres_0
293 sig(i,2) = -pres_0
294 sig(i,3) = -pres_0
295 ENDDO
296
297 IF(is_iform12)THEN
298
299 rho0_mix = vf(1)*rho0_(1)+vf(2)*rho0_(2)+vf(3)*rho0_(3)+vf(4)*rho0_(4)
300 e0_mix = vf(1)*e0_(1)+vf(2)*e0_(2)+vf(3)*e0_(3)+vf(4)*e0_(4)
301 !Sound Speed
302 k_mix = zero
303 if (rho0_(1) > zero .and. ssp_(1) > zero) k_mix = k_mix + vf(1) / ( rho0_(1) * ssp_(1)**2 )
304 if (rho0_(2) > zero .and. ssp_(2) > zero) k_mix = k_mix + vf(2) / ( rho0_(2) * ssp_(2)**2 )
305 if (rho0_(3) > zero .and. ssp_(3) > zero) k_mix = k_mix + vf(3) / ( rho0_(3) * ssp_(3)**2 )
306 if (rho0_(4) > zero .and. ssp_(4) > zero) k_mix = k_mix + vf(4) / ( rho0_(4) * ssp_(4)**2 )
307 If(k_mix > zero)k_mix=one/k_mix
308 ssp_mix=zero
309 if(gbuf%rho(i) > zero) ssp_mix = sqrt(k_mix / rho0_mix)
310 pm(27,mat(1)) = ssp_mix
311
312 DO i=1,nel
313 gbuf%RHO(i) = rho0_mix
314 gbuf%EINT(i) = e0_mix
315 !Temperature
316 v1=gbuf%VOL(i)*vf(1)
317 v2=gbuf%VOL(i)*vf(2)
318 v3=gbuf%VOL(i)*vf(3)
319 v4=gbuf%VOL(i)*vf(4)
320 denom = (sph1*v1+sph2*v2+sph3*v3+sph4*v4)
321 IF(denom /= zero)THEN
322 gbuf%TEMP(i) = (sph1*v1*t0_(1)+sph2*v2*t0_(2)+sph3*v3*t0_(3)+sph4*v4*t0_(4)) / denom
323 ENDIF
324 ENDDO
325 ENDIF
326
327 !---------------------------------!
328 ! IFORM 6 !
329 ! Initialize initial def param !
330 !---------------------------------!
331 IF(iform == 6) THEN
332 nv46 = 4
333 numel=numelq+numeltg
334 IF(n2d==0)THEN
335 nv46 = 6
336 numel=numels
337 ENDIF
338 CALL nrf51ini (
339 . ipm , pm , x , nix , ix,
340 . ale_connectivity , bufmat, uparam, rho0 ,
341 . uvar , nuvar , nel , gbuf%RHO, numel
342 . )
343 ENDIF
344C---
345 RETURN
346 END SUBROUTINE m51init
m51init
subroutine m51init(ipm, detonators, pm, tb, nuvar, uvar, uparam, x, mat, iparg, iform, ix, nix, ale_connectivity, bufmat, rho0, gbuf, nel, sig, mat_param, npf, tf)
Definition m51init.F:40
m5in2
subroutine m5in2(pm, mat, m151_id, detonators, tb, x, ix, nix)
Definition m5in2.F:40
m5in3
subroutine m5in3(pm, mat, m151_id, detonators, tb, iparg, x, ix, nix)
Definition m5in3.F:39
max
#define max(a, b)
Definition macros.h:21
ale_connectivity_mod
Definition ale_connectivity_mod.F:225
detonators_mod
Definition detonators_mod.F:80
eosmain_mod
Definition eosmain.F:35
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:102
nrf51ini
subroutine nrf51ini(ipm, pm, x, nix, ix, ale_connectivity, bufmat, uparam, rho0, uvar, nuvar, nel, rho, numel)
Definition nrf51ini.F:36
starter
program starter
Definition starter.F:39
ale_connectivity_mod::t_ale_connectivity
Definition ale_connectivity_mod.F:256
detonators_mod::detonators_struct_
Definition detonators_mod.F:155