40
41
42
43 USE elbufdef_mod
46 USE multimat_param_mod , ONLY : m51_n0phas, m51_nvphas
47 USE matparam_def_mod , ONLY : matparam_struct_
49
50
51
52#include "implicit_f.inc"
53
54
55
56#include "com01_c.inc"
57#include "com04_c.inc"
58#include "param_c.inc"
59
60
61
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) ::
70 TYPE (MATPARAM_STRUCT_), DIMENSION(NUMMAT), INTENT(INOUT) :: MAT_PARAM
71 INTEGER,INTENT(IN) :: NPF(*)
73
74
75
76 INTEGER :: I,, 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
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
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
94 LOGICAL :: IS_IFORM12
95 integer, parameter :: nvartmp = 3
96 integer, parameter :: nvareos = 6
97 my_real,
DIMENSION(1,nvareos) :: vareos
98 integer :: vartmp(1, nvartmp)
99
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
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
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
130 iparg(16) = 1
131 iparg(63) = 1
132 iparg(64) = 0
133 ENDIF
134
135
136
137
138 IF(iform == 2 .OR. iform == 3 .OR. iform == 4 .OR. iform==THEN
139 iparg(15) = 0
140 iparg(16) = 0
141 iparg(63) = 1
142 iparg(64) = 1
143 ENDIF
144
145
146
147
148 IF(iflg == 1 .AND. iparg(64) == 0)THEN
149 nph = 1
150 IF(n2d == 0)THEN
151 CALL m5in3 (pm,mat,0,detonators,tb,iparg,x
152 ELSE
153 CALL m5in2 (pm,mat,0,detonators,tb,x,ix,nix)
154 ENDIF
155 ENDIF
156
157
158
159
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
176
177 IF(.NOT. is_iform12)THEN
178
179
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))
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
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
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
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
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
244
245 ENDIF
246
247
248
249
250
251
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)
265 uvar(i,9+kk) = rho0_(imat)
266 uvar(i,10+kk) = zero
267 uvar(i,11+kk) = vold
268 uvar(i,12+kk) = rho0_(imat)
269 uvar(i,14+kk) = ssp_(imat)
270 uvar(i,15+kk) = zero
271 uvar(i,16+kk) = t0_(imat)
272 uvar(i,17+kk) = zero
273 uvar(i,18+kk) = p0_(imat)
274 uvar(i,19+kk) = zero
275 uvar(i,20+kk) = rho0_(imat)
276 uvar(i,21+kk) = e0_(imat)
277 uvar(i,22+kk) = zero
278 uvar(i,23+kk) = vf(imat)
279
280 ENDDO
281 ENDDO
282
283 imat = 4
284 kk = m51_n0phas+(imat-1)*m51_nvphas
285 uvar(1:nel,15+kk) = uvar(1:nel,1)
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
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
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)
323 ENDIF
324 ENDDO
325 ENDIF
326
327
328
329
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
339 . ipm , pm , x
340 . ale_connectivity , bufmat, uparam, rho0 ,
341 . uvar , nuvar , nel , gbuf%RHO, numel
342 . )
343 ENDIF
344
345 RETURN
subroutine m5in2(pm, mat, m151_id, detonators, tb, x, ix, nix)
subroutine m5in3(pm, mat, m151_id, detonators, tb, iparg, x, ix, nix)
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)
subroutine nrf51ini(ipm, pm, x, nix, ix, ale_connectivity, bufmat, uparam, rho0, uvar, nuvar, nel, rho, numel)