OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps97.F File Reference
#include "implicit_f.inc"
#include "comlock.inc"
#include "param_c.inc"
#include "com01_c.inc"
#include "com08_c.inc"
#include "mvsiz_p.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine sigeps97 (nel, nuparam, nuvar, tburn, time, uparam, bfrac, rho0, rho, eint, deltax, epspxx, epspyy, epspzz, signxx, signyy, signzz, signxy, signyz, signzx, sigvxx, sigvyy, sigvzz, sigvxy, sigvyz, sigvzx, soundsp, viscmax, uvar, off, geo, pid, ilay, ng, elbuf_tab, voln, qnew, qold, dpde)

Function/Subroutine Documentation

◆ sigeps97()

subroutine sigeps97 ( integer nel,
integer nuparam,
integer nuvar,
tburn,
time,
uparam,
bfrac,
rho0,
rho,
eint,
deltax,
epspxx,
epspyy,
epspzz,
signxx,
signyy,
signzz,
signxy,
signyz,
signzx,
sigvxx,
sigvyy,
sigvzz,
sigvxy,
sigvyz,
sigvzx,
soundsp,
viscmax,
uvar,
off,
geo,
integer, dimension(*) pid,
integer ilay,
integer ng,
type (elbuf_struct_), dimension(ngroup), target elbuf_tab,
voln,
qnew,
qold,
dimension(nel), intent(inout) dpde )

Definition at line 32 of file sigeps97.F.

43C-----------------------------------------------
44C M o d u l e s
45C-----------------------------------------------
46 USE elbufdef_mod
47 USE i22bufbric_mod
48 USE i22tri_mod
49C-----------------------------------------------
50C I m p l i c i t T y p e s
51C-----------------------------------------------
52#include "implicit_f.inc"
53#include "comlock.inc"
54C-----------------------------------------------
55C C o m m o n B l o c k s
56C-----------------------------------------------
57#include "param_c.inc"
58#include "com01_c.inc"
59#include "com08_c.inc"
60#include "mvsiz_p.inc"
61C-----------------------------------------------
62C I N P U T A r g u m e n t s
63C-----------------------------------------------
64 INTEGER NEL, NUPARAM, NUVAR, PID(*), ILAY, NG
65 my_real
66 . time ,uparam(nuparam),
67 . rho(nel) ,rho0(nel) ,
68 . eint(nel) ,qnew(nel) ,
69 . qold(nel) ,
70 . epspxx(nel) ,epspyy(nel),epspzz(nel),
71 . geo(npropg,*),ssp ,
72 . voln(*)
73 TYPE (ELBUF_STRUCT_), TARGET, DIMENSION(NGROUP) :: ELBUF_TAB
74C-----------------------------------------------
75C O U T P U T A r g u m e n t s
76C-----------------------------------------------
77 my_real
78 . signxx(nel),signyy(nel),signzz(nel),
79 . signxy(nel),signyz(nel),signzx(nel),
80 . sigvxx(nel),sigvyy(nel),sigvzz(nel),
81 . sigvxy(nel),sigvyz(nel),sigvzx(nel),
82 . soundsp(nel),viscmax(nel)
83 my_real,INTENT(INOUT) :: dpde(nel)
84C-----------------------------------------------
85C I N P U T O U T P U T A r g u m e n t s
86C-----------------------------------------------
87 my_real uvar(nel,nuvar), off(nel), tburn(nel), bfrac(nel), deltax(nel)
88C-----------------------------------------------
89C VARIABLES FOR FUNCTION INTERPOLATION
90C-----------------------------------------------
91C-----------------------------------------------
92C L o c a l V a r i a b l e s
93C-----------------------------------------------
94c my_real
95c . B1, B2, R1, R2, W1, D, PCJ, E0, C1, VCJ,
96c . EADD, TBEGIN, TEND,BHE,
97c . PSH,REACTION_RATE,REACTION_RATE2,A_MIL,M_MIL,N_MIL,ALPHA_UNIT,
98c . DPDMU,
99c . R1V, R2V, ER1M, ER2M, MUP1, MU, AA,BB, P0, P, VOLD, POLD,
100c . WDR1V,WDR2V,DR1V,DR2V,ER1V,ER2V, PNEW,VOLO,FACM,ESPE,DVOL,
101c . W1DF,EINC,QA,QB,QAL,QBL,DD
102c INTEGER :: IBID, IBFRAC, QOPT, I
103c my_real :: XL, MASS, DF, TB
104
105 my_real
106 . d, pcj, e0, p0, vcj,c,psh,
107 . a(5),r(5),al(5),bl(5),rl(5),
108 . pnew,espe,dvol,
109 . qa,qb,qal,qbl,dd,bhe,
110 . lambda1,lambda2,lambda3,lambda4,lambda5,
111 . lambda,
112 . dldv1,dldv2,dldv3,dldv4,dldv5,
113 . dldv,
114 . erlv1,erlv2,erlv3,erlv4,erlv5,
115 . p1,p2,p3,p4,p5,
116 . rhoc2,
117 . rv1,rv2,rv3,rv4,rv5,
118 . rhoc2_1,rhoc2_2,rhoc2_3,rhoc2_4,rhoc2_5,
119 . dv, vv, vv_,
120 . ww, denom, p, v0
121
122 INTEGER :: IBFRAC, I
123 my_real :: xl, df, tb
124
125 TYPE(G_BUFEL_) ,POINTER :: GBUF
126 TYPE(L_BUFEL_) ,POINTER :: LBUF
127 TYPE(BUF_LAY_) ,POINTER :: BUFLY
128C-----------------------------------------------
129C S o u r c e L i n e s
130C-----------------------------------------------
131 gbuf => elbuf_tab(ng)%GBUF
132 lbuf => elbuf_tab(ng)%BUFLY(ilay)%LBUF(1,1,1)
133 bufly => elbuf_tab(ng)%BUFLY(ilay)
134
135
136C-----------------------------------------------
137C S t a n d a r d J W L E O S
138C-----------------------------------------------
139 !IN CASE OF SWITCHING LAW5 LATER FROM LECMAT/MMAIN TO LECMUSER/MULAW
140
141c E0 = UPARAM(1)
142c P0 = UPARAM(2)
143c B1 = UPARAM(4)
144c B2 = UPARAM(5)
145c R1 = UPARAM(6)
146c R2 = UPARAM(7)
147c W1 = UPARAM(8)
148c D = UPARAM(9)
149c PCJ = UPARAM(10)
150c BHE = UPARAM(11)
151c IBFRAC = UPARAM(12)
152c QOPT = UPARAM(13)
153c EADD = UPARAM(14)
154c TBEGIN = UPARAM(15)
155c TEND = UPARAM(16)
156c REACTION_RATE = UPARAM(17)
157c A_MIL = UPARAM(18)
158c M_MIL = UPARAM(19)
159c N_MIL = UPARAM(20)
160c REACTION_RATE2 = UPARAM(21)
161c ALPHA_UNIT = UPARAM(22)
162c PSH = UPARAM(23)
163c
164c IF(DT1==ZERO)THEN
165c DO I=1,NEL
166c EINT(I) = E0*VOLN(I)
167c ENDDO
168c UVAR(1:NEL,4) = VOLN(1:NEL) !VOLD
169c UVAR(1:NEL,5) = P0 !POLD
170c ENDIF
171c DO I=1,NEL
172c !--------------------------------!
173c ! Calculation of BFRAC in [0,1] !
174c !--------------------------------!
175c XL = DELTAX(I) !VOLN(I)**THIRD
176c IF(BFRAC(I) < ONE) THEN
177c TB = - TBURN(I)
178c BFRAC(I) = ZERO
179c IF(IBFRAC/=1 .AND. TIME > TB) BFRAC(I) = D*(TIME-TB)*TWO_THIRD/XL
180c IF(IBFRAC/=2) BFRAC(I) = MAX( BFRAC(I) , BHE * (ONE - RHO0(I)/RHO(I)) )
181c IF(BFRAC(I) < EM04) THEN
182c BFRAC(I) = ZERO
183c ELSEIF(BFRAC(I) > ONE) THEN
184c BFRAC(I) = ONE
185c ENDIF
186c ENDIF
187c !--------------------------------!
188c ! EOS SOLVING !
189c !--------------------------------!
190c DF = RHO0(I)/RHO(I)
191c R1V = B1*W1/(R1*DF)
192c R2V = B2*W1/(R2*DF)
193c WDR1V = B1-R1V
194c WDR2V = B2-R2V
195c DR1V = W1*EINT(I)/MAX(EM20,VOLN(I)) !w*Eint/V = w*E/v where v=V/V0
196c ER1V = EXP(-R1*DF)
197c ER2V = EXP(-R2*DF)
198c P = P0 - PSH + WDR1V*ER1V+WDR2V*ER2V+DR1V
199c P = MAX(ZERO - PSH, P)
200c SSP = B1*ER1V*( (-W1/DF/R1) + R1*DF - W1)
201c . + B2*ER2V*( (-W1/DF/R2) + R2*DF - W1)
202c . + DR1V + (P + PSH)*W1
203c SSP = SSP * DF
204c SSP = SQRT(ABS(SSP)/RHO0(I))
205c SSP = MAX(SSP,D*(ONE-BFRAC(I)))
206c QA = GEO(14,PID(I))
207c QB = GEO(15,PID(I))
208c DD = -EPSPXX(I)-EPSPYY(I)-EPSPZZ(I)
209c QAL = QA*XL
210c QAL = QAL*QAL
211c QBL = QB*XL
212c VISCMAX(I) = RHO(I)*(QAL*MAX(ZERO,DD) + QBL*SSP)
213c QNEW(I) = VISCMAX(I)*MAX(ZERO,DD)
214c viscmax(I) = zero
215c !QNEW(I) = zero !QOLD(I)
216c VOLD = UVAR(I,4)
217c POLD = -UVAR(I,5)
218c DVOL = HALF*(VOLN(I)-VOLD)
219c EINC = DVOL*(POLD-PSH-QOLD(I)-QNEW(I))
220c !EINT(I) = EINT(I)+EINC
221c QOLD(I) = QNEW(I)
222c VOLO = VOLN(I)/DF
223c ESPE = (EINT(I)+EINC)/MAX(EM20,VOLO)
224c W1DF = BFRAC(I)*W1/DF
225c FACM = BFRAC(I)*(WDR1V*ER1V+WDR2V*ER2V)
226c PNEW = P0 - PSH + (FACM+ESPE*W1DF)/(ONE +W1DF*DVOL/MAX(EM20,VOLO))
227c PNEW = MAX(ZERO - PSH, PNEW)*OFF(I)
228c EINC = EINC-(PNEW + PSH)*DVOL
229c ! EINT(I) = (EINT(I)-(PNEW + PSH)*DVOL)!/MAX(EM20,VOLN(I))
230c UVAR(I,4) = VOLN(I)
231c UVAR(I,5) = PNEW
232c !--------------------------------!
233c ! Returning values !
234c !--------------------------------!
235c SIGNXX(I) = -PNEW
236c SIGNYY(I) = -PNEW
237c SIGNZZ(I) = -PNEW
238c SIGNXY(I) = ZERO
239c SIGNYZ(I) = ZERO
240c SIGNZX(I) = ZERO
241c SIGVXX(I) = ZERO
242c SIGVYY(I) = ZERO
243c SIGVZZ(I) = ZERO
244c SIGVXY(I) = ZERO
245c SIGVYZ(I) = ZERO
246c SIGVZX(I) = ZERO
247c
248c SOUNDSP(I) = SSP
249c ENDDO! next I
250
251
252C-----------------------------------------------
253C J W L - B A K E R E O S
254C-----------------------------------------------
255
256 p0 = uparam(01)
257 psh = uparam(02)
258 ibfrac = nint(uparam(03))
259 d = uparam(04)
260 pcj = uparam(05)
261 e0 = uparam(06)
262 ww = uparam(07)
263 c = uparam(08)
264 a(1:5) = uparam(09:13)
265 r(1:5) = uparam(14:18)
266 al(1:5)= uparam(19:23)
267 bl(1:5)= uparam(24:28)
268 rl(1:5)= uparam(29:33)
269 bhe = uparam(34)
270 vcj = uparam(35)
271
272
273 IF(dt1==zero)THEN
274 DO i=1,nel
275 eint(i) = e0*voln(i)
276 ENDDO
277 uvar(1:nel,4) = voln(1:nel) !VOLD
278 uvar(1:nel,5) = p0 !POLD
279 ENDIF
280
281 DO i=1,nel
282 !--------------------------------!
283 ! Calculation of BFRAC in [0,1] !
284 !--------------------------------!
285 xl = deltax(i) !VOLN(I)**THIRD
286 IF(bfrac(i) < one) THEN
287 tb = - tburn(i)
288 bfrac(i) = zero
289 IF(ibfrac/=1 .AND. time > tb) bfrac(i) = d*(time-tb)*two_third/xl
290 IF(ibfrac/=2) bfrac(i) = max( bfrac(i) , bhe * (one - rho0(i)/rho(i)) )
291 IF(bfrac(i) < em04) THEN
292 bfrac(i) = zero
293 ELSEIF(bfrac(i) > one) THEN
294 bfrac(i) = one
295 ENDIF
296 ENDIF
297 !--------------------------------!
298 ! EOS SOLVING !
299 !--------------------------------!
300 df = rho0(i)/rho(i)
301 v0 = rho(i)*voln(i) / rho0(i)
302 espe = eint(i)/max(em20,v0)
303
304 ! REMINDER
305 !UVAR(I,4) = VOLN(I)
306 !UVAR(I,5) = PNEW
307
308
309 vv = voln(i) / v0 ! current relative volume
310 vv_ = uvar(i,4) / v0 ! previous relative volume
311 dv = vv-vv_
312
313 erlv1 = exp(-rl(1)*vv)
314 erlv2 = exp(-rl(2)*vv)
315 erlv3 = exp(-rl(3)*vv)
316 erlv4 = exp(-rl(4)*vv)
317 erlv5 = exp(-rl(5)*vv)
318
319 lambda1 = (al(1)*vv+bl(1))*erlv1
320 lambda2 = (al(2)*vv+bl(2))*erlv2
321 lambda3 = (al(3)*vv+bl(3))*erlv3
322 lambda4 = (al(4)*vv+bl(4))*erlv4
323 lambda5 = (al(5)*vv+bl(5))*erlv5
324
325 lambda = lambda1 + lambda2 + lambda3 + lambda4 + lambda5 + ww
326
327 dpde(i) = lambda/df
328
329 dldv1 = al(1)*erlv1-(al(1)*vv+bl(1))*rl(1)*erlv1
330 dldv2 = al(2)*erlv2-(al(2)*vv+bl(2))*rl(2)*erlv2
331 dldv3 = al(3)*erlv3-(al(3)*vv+bl(3))*rl(3)*erlv3
332 dldv4 = al(4)*erlv4-(al(4)*vv+bl(4))*rl(4)*erlv4
333 dldv5 = al(5)*erlv5-(al(5)*vv+bl(5))*rl(5)*erlv5
334
335 dldv = dldv1 + dldv2 + dldv3 + dldv4 + dldv5
336
337 rv1 = r(1)*vv
338 rv2 = r(2)*vv
339 rv3 = r(3)*vv
340 rv4 = r(4)*vv
341 rv5 = r(5)*vv
342
343 p1 = a(1)*(one-lambda/rv1)*exp(-rv1)
344 p2 = a(2)*(one-lambda/rv2)*exp(-rv2)
345 p3 = a(3)*(one-lambda/rv3)*exp(-rv3)
346 p4 = a(4)*(one-lambda/rv4)*exp(-rv4)
347 p5 = a(5)*(one-lambda/rv5)*exp(-rv5)
348
349 p = p1+p2+p3+p4+p5 + lambda*espe/vv + c*(one-lambda/ww)*exp((-ww-one)*log(vv))
350
351 rhoc2_1 = a(1)*( (vv*dldv-lambda)/r(1) + r(1)*vv*vv - lambda*vv )*exp(-rv1)
352 rhoc2_2 = a(2)*( (vv*dldv-lambda)/r(2) + r(2)*vv*vv - lambda*vv )*exp(-rv2)
353 rhoc2_3 = a(3)*( (vv*dldv-lambda)/r(3) + r(3)*vv*vv - lambda*vv )*exp(-rv3)
354 rhoc2_4 = a(4)*( (vv*dldv-lambda)/r(4) + r(4)*vv*vv - lambda*vv )*exp(-rv4)
355 rhoc2_5 = a(5)*( (vv*dldv-lambda)/r(5) + r(5)*vv*vv - lambda*vv )*exp(-rv5)
356
357 !EINC = -HALF*(UVAR(I,5)+P+PSH+PSH)*Dv
358
359 rhoc2 = rhoc2_1 + rhoc2_2 + rhoc2_3 + rhoc2_4 + rhoc2_5
360 rhoc2 = rhoc2 + c*((ww+one)*(one-lambda/ww)+vv*dldv/ww)*exp(-ww*log(vv))
361 rhoc2 = rhoc2 + (espe)*lambda + lambda*vv*(p+psh) - (espe)*vv*dldv
362
363 ssp = sqrt(max(rhoc2/rho0(i),em20))
364 qa = geo(14,pid(i))
365 qb = geo(15,pid(i))
366 dd = -epspxx(i)-epspyy(i)-epspzz(i)
367 qal = qa*xl
368 qal = qal*qal
369 qbl = qb*xl
370 viscmax(i) = rho(i)*(qal*max(zero,dd) + qbl*ssp)
371 qnew(i) = viscmax(i)*max(zero,dd)
372 !viscmax(I) = zero
373
374 denom = (one+half*dv*lambda/vv)
375 denom = one/denom
376
377 !CURRENT PRESSURE
378 pnew = p - lambda/vv*half*(uvar(i,5)+psh)*dv
379 pnew = denom * pnew
380 pnew = (one-bfrac(i))*p0 + bfrac(i)*pnew
381 pnew = max(-psh, pnew-psh)*off(i)
382
383 !CURRENT ENERGY
384 !pold = uvar(i,5)
385 dvol = half*(voln(i)-uvar(i,4))
386 !EINC = DVOL*(-PNEW-PSH-POLD-PSH-QOLD(I)-QNEW(I))
387 eint(i) = eint(i) - (psh+psh)*dvol
388
389 !BACKUP FOR NEXT CYCLE
390 qold(i) = qnew(i)
391 uvar(i,4) = voln(i)
392 uvar(i,5) = pnew
393
394 !--------------------------------!
395 ! Returning values !
396 !--------------------------------!
397
398 signxx(i) = -pnew
399 signyy(i) = -pnew
400 signzz(i) = -pnew
401 signxy(i) = zero
402 signyz(i) = zero
403 signzx(i) = zero
404
405 sigvxx(i) = zero
406 sigvyy(i) = zero
407 sigvzz(i) = zero
408 sigvxy(i) = zero
409 sigvyz(i) = zero
410 sigvzx(i) = zero
411
412 soundsp(i) = ssp
413
414 enddo! next I
415
416
417C-----------------------------------------------
418 RETURN
419C-----------------------------------------------
my_real
#define my_real
Definition cppsort.cpp:32
max
#define max(a, b)
Definition macros.h:21
i22bufbric_mod
Definition cut-cell-search_mod.F:95
i22tri_mod
Definition cut-cell-search_mod.F:212