OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps75.F File Reference
#include "implicit_f.inc"
#include "com04_c.inc"
#include "param_c.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine sigeps75 (nel, nuparam, nuvar, nfunc, ifunc, npf, tf, time, timestep, uparam, rho0, rho, vnew, eint, muold, epspxx, epspyy, epspzz, epspxy, epspyz, epspzx, depsxx, depsyy, depszz, depsxy, depsyz, depszx, epsxx, epsyy, epszz, epsxy, epsyz, epszx, sigoxx, sigoyy, sigozz, sigoxy, sigoyz, sigozx, signxx, signyy, signzz, signxy, signyz, signzx, sigvxx, sigvyy, sigvzz, sigvxy, sigvyz, sigvzx, soundsp, viscmax, uvar, off, dvol, vol0, pm, ipm, mat, psh, bufmat, vareos, nvareos, mat_param, nvartmp_eos, vartmp_eos)

Function/Subroutine Documentation

◆ sigeps75()

subroutine sigeps75 ( integer nel,
integer nuparam,
integer nuvar,
integer nfunc,
integer, dimension(nfunc) ifunc,
integer, dimension(*) npf,
tf,
time,
timestep,
uparam,
rho0,
rho,
vnew,
eint,
muold,
epspxx,
epspyy,
epspzz,
epspxy,
epspyz,
epspzx,
depsxx,
depsyy,
depszz,
depsxy,
depsyz,
depszx,
epsxx,
epsyy,
epszz,
epsxy,
epsyz,
epszx,
sigoxx,
sigoyy,
sigozz,
sigoxy,
sigoyz,
sigozx,
signxx,
signyy,
signzz,
signxy,
signyz,
signzx,
sigvxx,
sigvyy,
sigvzz,
sigvxy,
sigvyz,
sigvzx,
soundsp,
viscmax,
uvar,
off,
dvol,
vol0,
pm,
integer, dimension(npropmi,numgeo) ipm,
integer, dimension(nel) mat,
psh,
bufmat,
dimension(nvareos*nel), intent(inout) vareos,
integer, intent(in) nvareos,
type(matparam_struct_), dimension(nummat), intent(in) mat_param,
integer, intent(in) nvartmp_eos,
integer, dimension(nel,nvartmp_eos), intent(inout) vartmp_eos )

Definition at line 37 of file sigeps75.F.

50C
51C---------+---------+---+---+--------------------------------------------
52C VAR | SIZE |TYP| RW| DEFINITION
53C---------+---------+---+---+--------------------------------------------
54C NEL | 1 | I | R | SIZE OF THE ELEMENT GROUP NEL
55C NUPARAM | 1 | I | R | SIZE OF THE USER PARAMETER ARRAY
56C NUVAR | 1 | I | R | NUMBER OF USER ELEMENT VARIABLES
57C---------+---------+---+---+--------------------------------------------
58C MFUNC | 1 | I | R | NUMBER FUNCTION USED FOR THIS USER LAW not used
59C KFUNC | NFUNC | I | R | FUNCTION INDEX not used
60C NPF | * | I | R | FUNCTION ARRAY
61C TF | * | F | R | FUNCTION ARRAY
62C---------+---------+---+---+--------------------------------------------
63C TIME | 1 | F | R | CURRENT TIME
64C TIMESTEP| 1 | F | R | CURRENT TIME STEP
65C UPARAM | NUPARAM | F | R | USER MATERIAL PARAMETER ARRAY
66C RHO0 | NEL | F | R | INITIAL DENSITY
67C RHO | NEL | F | R | DENSITY
68C VNEW | NEL | F | R | VOLUME
69C VOL0 | NEL | F | R | VOLUME INITIAL (LAGRANGE)
70C | | F | R | VOLUME (ALE)
71C EINT | NEL | F | R | TOTAL INTERNAL ENERGY
72C EPSPXX | NEL | F | R | STRAIN RATE XX
73C EPSPYY | NEL | F | R | STRAIN RATE YY
74C ... | | | |
75C DEPSXX | NEL | F | R | STRAIN INCREMENT XX
76C DEPSYY | NEL | F | R | STRAIN INCREMENT YY
77C ... | | | |
78C EPSXX | NEL | F | R | STRAIN XX
79C EPSYY | NEL | F | R | STRAIN YY
80C ... | | | |
81C SIGOXX | NEL | F | R | OLD ELASTO PLASTIC STRESS XX
82C SIGOYY | NEL | F | R | OLD ELASTO PLASTIC STRESS YY
83C ... | | | |
84C---------+---------+---+---+--------------------------------------------
85C SIGNXX | NEL | F | W | NEW ELASTO PLASTIC STRESS XX
86C SIGNYY | NEL | F | W | NEW ELASTO PLASTIC STRESS YY
87C ... | | | |
88C SIGVXX | NEL | F | W | VISCOUS STRESS XX
89C SIGVYY | NEL | F | W | VISCOUS STRESS YY
90C ... | | | |
91C SOUNDSP | NEL | F | W | SOUND SPEED (NEEDED FOR TIME STEP)
92C VISCMAX | NEL | F | W | MAXIMUM DAMPING MODULUS(NEEDED FOR TIME STEP)
93C---------+---------+---+---+--------------------------------------------
94C UVAR |NEL*NUVAR| F |R/W| USER ELEMENT VARIABLE ARRAY
95C OFF | NEL | F |R/W| DELETED ELEMENT FLAG (=1. ON, =0. OFF)
96C---------+---------+---+---+--------------------------------------------
97C-----------------------------------------------
98C M o d u l e s
99C-----------------------------------------------
100 USE message_mod
101 USE matparam_def_mod
102 USE eosmain_mod , ONLY : eosmain
103C-----------------------------------------------
104C I M P L I C I T T Y P E S
105C-----------------------------------------------
106#include "implicit_f.inc"
107C-----------------------------------------------
108C C O M M O N
109C-----------------------------------------------
110#include "com04_c.inc"
111#include "param_c.inc"
112C----------------------------------------------------------------
113C I N P U T A R G U M E N T S
114C----------------------------------------------------------------
115 INTEGER NEL, NUPARAM, NUVAR, IPM(NPROPMI,NUMGEO), MAT(NEL)
116 INTEGER,INTENT(IN) ::NVAREOS
117 my_real
118 . time , timestep , uparam(nuparam),
119 . rho(nel), rho0(nel), vnew(nel), eint(nel),
120 . epspxx(nel), epspyy(nel), epspzz(nel),
121 . epspxy(nel), epspyz(nel), epspzx(nel),
122 . depsxx(nel), depsyy(nel), depszz(nel),
123 . depsxy(nel), depsyz(nel), depszx(nel),
124 . epsxx(nel), epsyy(nel), epszz(nel),
125 . epsxy(nel), epsyz(nel), epszx(nel),
126 . sigoxx(nel), sigoyy(nel), sigozz(nel),
127 . sigoxy(nel), sigoyz(nel), sigozx(nel),
128 . dvol(nel) , vol0(nel) , pm(npropm,nummat),
129 . psh(nel) , bufmat(*) , muold(nel)
130
131 my_real,INTENT(INOUT) :: vareos(nvareos*nel)
132 TYPE(MATPARAM_STRUCT_) ,DIMENSION(NUMMAT), INTENT(IN) :: MAT_PARAM
133 INTEGER,INTENT(IN) :: NVARTMP_EOS
134 INTEGER,INTENT(INOUT) :: VARTMP_EOS(NEL,NVARTMP_EOS)
135C----------------------------------------------------------------
136C O U T P U T A R G U M E N T S
137C----------------------------------------------------------------
138 my_real
139 . signxx(nel), signyy(nel), signzz(nel),
140 . signxy(nel), signyz(nel), signzx(nel),
141 . sigvxx(nel), sigvyy(nel), sigvzz(nel),
142 . sigvxy(nel), sigvyz(nel), sigvzx(nel),
143 . soundsp(nel), viscmax(nel)
144C----------------------------------------------------------------
145C I N P U T O U T P U T A R G U M E N T S
146C----------------------------------------------------------------
147 my_real uvar(nel,nuvar), off(nel)
148C----------------------------------------------------------------
149C VARIABLES FOR FUNCTION INTERPOLATION
150C----------------------------------------------------------------
151 INTEGER NPF(*), NFUNC, IFUNC(NFUNC)
152 my_real finter,fintte,tf(*),fint2v
153 EXTERNAL finter,fintte
154C----------------------------------------------------------------
155C L o c a l V a r i a b l e s
156C----------------------------------------------------------------
157 INTEGER I, ITER, MATS(NEL),IPLA(NEL), MXS, IFLAG1,
158 . IFLAG2, ITEMAX, IMAX, IPLAS, EOSTYP
159 my_real :: bid(1) !bfrac not used in this context
160 my_real
161 . dfs(nel), mus(nel), mu2s(nel), dpdm(nel),
162 . theta(nel), dpde(nel)
163 my_real
164 . mu(nel)
165 my_real
166 . ps,pe,alphae,alphap,en,beta,eta,xi,unsurn,
167 . bulk0,c0,dalpdpe,aa,hh,hh2,g,g2,dav,
168 . bulk1, crit, tol, pnew1, dfda, dgda,
169 . gs, dgdalpha, bulk
170 my_real
171 . rhos(nel),rhos0(nel),espes(nel),
172 . vnews(nel),dvols(nel),pnew(nel),
173 . alphayld(nel),pyld(nel),alpha0(nel),
174 . alphaold(nel),dpdael(nel),eintn(nel),pold(nel)
175 my_real
176 . sigold(nel,6), pnewg(nel)
177 DOUBLE PRECISION ALPHA1(NEL),ALPHA(NEL),DALPHA
178C---+----1----+----2----+----3----+----4----+----5----+----6----+----7--
179 bulk=uparam(1)
180 pe = uparam(2)
181 ps = uparam(3)
182 en = uparam(4)
183 tol= uparam(5)
184 mxs =nint(uparam(6))
185 iflag1=nint(uparam(7))
186 iflag2=nint(uparam(8))
187 itemax=nint(uparam(17))
188 DO i=1,nel
189 mats(i) = mxs
190 rhos0(i)= pm(89,mxs)
191 ENDDO
192 bulk0= pm(32,mxs)
193 g =uparam(10)
194 alphae=uparam(11)
195 gs = pm(22,mxs)
196 dgdalpha= (g-gs)/(alphae-one)
197 alphap=uparam(12)
198 aa =uparam(13)
199 unsurn=uparam(14)
200 xi =uparam(15)
201 eta =uparam(16)
202 eostyp = mat_param(mxs)%IEOS ! LAW 75 uses /EOS defined in another law !!!
203C
204 DO i=1,nel
205 alphaold(i)= uvar(i,3)
206 alphayld(i)= uvar(i,4)
207 pold(i)=-(sigoxx(i)+sigoyy(i)+sigozz(i))*third
208 pyld(i)=ps-xi*(alphayld(i)-one)**unsurn
209 ENDDO
210 DO i=1,nel
211 mu(i) =rho(i)/rho0(i)-one
212 ENDDO
213 DO i=1,nel
214 sigold(i,1:3)=-pold(i)
215 sigold(i,4:6)=zero
216 ENDDO
217C----------------------------
218C FIRST GUESS LINEAR SOLUTION
219C----------------------------
220 DO i=1,nel
221 IF(alphaold(i) == one) THEN
222C COMPACT MATERIAL
223 alpha0(i)=one
224 ipla(i)=2
225 ELSE
226C POROUS MATERIAL
227 hh =aa*(alphaold(i)-one)+one
228 hh2=hh*hh
229 dpdael(i)=bulk0/(alphaold(i)*(one-alphaold(i)/hh2))
230 bulk1=bulk0*(one+mu(i))/alphae
231 IF(iflag1 == 1) THEN
232 alpha0(i)=(pold(i)+bulk0-alphaold(i)*dpdael(i))/(bulk1-dpdael(i))
233 ELSEIF(iflag1 == 2) THEN
234 alpha0(i)= bulk0/(bulk1-pold(i))
235 ENDIF
236 alpha0(i)=max(alpha0(i),one)
237 ENDIF
238 ENDDO
239C
240 DO i=1,nel
241 alpha(i)=alpha0(i)
242 ENDDO
243C---------------------
244C SOLUTION ELASTIQUE
245C---------------------
246 iter=0
247 100 iter=iter+1
248C----------------
249C SOLID VARIABLES
250C----------------
251 DO i=1,nel
252 rhos(i) =rho(i)*alpha(i)
253 dvols(i)=dvol(i)-vnew(i)*(one-alphaold(i)/alpha(i))
254 dvols(i)=dvols(i)/alphaold(i)
255 vnews(i)=vnew(i)/alpha(i)
256 ENDDO
257C------------------
258C ENERGY ESTIMATION
259C------------------
260 DO i=1,nel
261 mus(i) = (one + mu(i))*alpha(i)/alphae - one
262 dfs(i) = one/(one+mus(i))
263 eintn(i)=eint(i)-half*pold(i)*dvols(i)
264 ENDDO
265C--------------
266C EOS FOR SOLID
267C--------------
268 CALL eosmain(0 ,nel ,eostyp ,pm ,off ,eintn ,
269 . rhos ,rhos0 ,mus ,mu2s ,espes ,
270 . dvols ,dfs ,vnews ,mats ,psh ,
271 . pnew ,dpdm ,dpde ,theta ,
272 . bufmat ,sigold ,muold ,75 ,
273 . npf ,tf ,vareos ,nvareos,mat_param(mxs),
274 . bid ,nvartmp_eos, vartmp_eos)
275 CALL eosmain(1 ,nel ,eostyp ,pm ,off ,eintn ,
276 . rhos ,rhos0 ,mus ,mu2s ,espes ,
277 . dvols ,dfs ,vnews ,mats ,psh ,
278 . pnew ,dpdm ,dpde ,theta ,
279 . bufmat ,sigold ,muold ,75 ,
280 . npf ,tf ,vareos ,nvareos,mat_param(mxs),
281 . bid ,nvartmp_eos, vartmp_eos)
282 IF(iflag1 == 2) THEN
283 DO i=1,nel
284 pnew(i)=pnew(i)/alpha(i)
285 ENDDO
286 ENDIF
287C------------
288C CONVERGENCE
289C------------
290 IF (iter <= itemax) THEN
291C-----------------
292C NEWTON ITERATION
293C-----------------
294 DO i=1,nel
295 IF(ipla(i)==2)cycle
296 ipla(i)=0
297 dfda=dpdm(i)*(one+mu(i))/alphae + dpde(i)*eint(i)/vol0(i)
298 dgda=dpdael(i)
299 pnewg(i)=pold(i)+dpdael(i)*(alpha(i)-alphaold(i))
300C
301 IF(iflag1 == 1) THEN
302 dalpha=(pnew(i)-pnewg(i))/(dfda-dgda)
303 ELSEIF(iflag1 == 2) THEN
304 dalpha=(pnew(i)-pnewg(i))/((dfda-pnew(i))/alpha(i)-dgda)
305 ENDIF
306 alpha1(i)=alpha(i)-dalpha
307 alpha1(i)=max(alpha1(i),one)
308 ENDDO
309
310 crit=-one
311 DO i=1,nel
312 IF(ipla(i)==2)cycle
313 beta=abs(alpha(i)-alpha1(i))
314 IF(beta > crit) THEN
315 crit=beta
316 imax=i
317 ENDIF
318 ENDDO
319C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----A----+----B----+----C----+
320C WRITE(6,'(5X,A,I4,8(A,E18.10))')'ELAS ITER=',ITER,'MUS=',MUS(IMAX),' ALPHA1=',ALPHA1(IMAX),' ALPHA=',ALPHA(IMAX),
321C . ' CRIT=',CRIT,' PNEW=',PNEW(IMAX),' PNEWG=',PNEWG(IMAX),' ALPHAOLD=',ALPHAOLD(IMAX)
322C
323 IF (crit > tol) THEN
324 DO i=1,nel
325 alpha(i)=alpha1(i)
326 ENDDO
327 GO TO 100
328 ENDIF
329 ELSE
330 CALL ancmsg(msgid=137,anmode=aninfo,i1=itemax)
331 CALL arret(2)
332 ENDIF
333C-----------------------
334C CHECK ELASTIC SOLUTION
335C-----------------------
336 iplas=0
337 DO i=1,nel
338 IF(ipla(i)==2)cycle
339 IF(pnew(i) > pyld(i)) THEN
340 ipla(i)=1
341 iplas=1
342 ENDIF
343 ENDDO
344 IF(iplas==0) GO TO 300
345C---------------------
346C SOLUTION PLASTIQUE
347C---------------------
348 DO i=1,nel
349 IF(ipla(i)==1) alpha(i)=alpha0(i)
350 ENDDO
351 iter=0
352 200 iter=iter+1
353C----------------
354C SOLID VARIABLES
355C----------------
356 DO i=1,nel
357 rhos(i) =rho(i)*alpha(i)
358 dvols(i)=dvol(i)-vnew(i)*(one-alphaold(i)/alpha(i))
359 dvols(i)=dvols(i)/alphaold(i)
360 vnews(i)=vnew(i)/alpha(i)
361 ENDDO
362C------------------
363C ENERGY ESTIMATION
364C------------------
365 DO i=1,nel
366 mus(i) = (one + mu(i))*alpha(i)/alphae - one
367 dfs(i) = one/(one+mus(i))
368 eintn(i)=eint(i)-half*pold(i)*dvols(i)
369 ENDDO
370C--------------
371C EOS FOR SOLID
372C--------------
373 CALL eosmain(0 ,nel ,eostyp ,pm ,off ,eintn ,
374 . rhos ,rhos0 ,mus ,mu2s ,espes ,
375 . dvols ,dfs ,vnews ,mats ,psh ,
376 . pnew ,dpdm ,dpde ,theta ,
377 . bufmat,sigold,muold ,75 ,
378 . npf ,tf ,vareos ,nvareos,mat_param(mxs),
379 . bid, nvartmp_eos, vartmp_eos)
380 CALL eosmain(1 ,nel ,eostyp ,pm ,off ,eintn ,
381 . rhos ,rhos0 ,mus ,mu2s ,espes ,
382 . dvols ,dfs ,vnews ,mats ,psh ,
383 . pnew ,dpdm ,dpde ,theta ,
384 . bufmat,sigold,muold ,75 ,
385 . npf ,tf ,vareos ,nvareos,mat_param(mxs),
386 . bid, nvartmp_eos, vartmp_eos)
387 IF(iflag1 == 2) THEN
388 DO i=1,nel
389 pnew(i)=pnew(i)/alpha(i)
390 ENDDO
391 ENDIF
392C------------
393C CONVERGENCE
394C------------
395 IF (iter <= itemax) THEN
396C-----------------
397C NEWTON ITERATION
398C-----------------
399 DO i=1,nel
400 IF(ipla(i)/=1)cycle
401 dfda=dpdm(i)*(one+mu(i))/alphae + dpde(i)*eint(i)/vol0(i)
402 dgda=-xi*unsurn*(max(alpha(i)-one,em20))**(unsurn-one)
403 pnewg(i)=ps-xi*(max(alpha(i)-one,em20))**unsurn
404C
405 IF(iflag1 == 1) THEN
406 dalpha=(pnew(i)-pnewg(i))/(dfda-dgda)
407 ELSEIF(iflag1 == 2) THEN
408 dalpha=(pnew(i)-pnewg(i))/((dfda-pnew(i))/alpha(i)-dgda)
409 ENDIF
410 alpha1(i)=alpha(i)-dalpha
411 alpha1(i)=max(alpha1(i),one)
412 ENDDO
413
414 crit=-one
415 DO i=1,nel
416 IF(ipla(i)/=1)cycle
417 beta=abs(alpha(i)-alpha1(i))
418 IF(beta > crit) THEN
419 crit=beta
420 imax=i
421 ENDIF
422 ENDDO
423C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----A----+----B----+----C----+
424C WRITE(6,'(5X,A,I4,8(A,E18.10))')'PLAS ITER=',ITER,'MUS=',MUS(IMAX),' ALPHA1=',ALPHA1(IMAX),' ALPHA=',ALPHA(IMAX),
425C . ' CRIT=',CRIT,' PNEW=',PNEW(IMAX),' PNEWG=',PNEWG(IMAX),' ALPHAOLD=',ALPHAOLD(IMAX)
426C
427 IF (crit > tol) THEN
428 DO i=1,nel
429 alpha(i)=alpha1(i)
430 ENDDO
431 GO TO 200
432 ENDIF
433 ELSE
434 CALL ancmsg(msgid=137,anmode=aninfo,i1=itemax)
435 CALL arret(2)
436 ENDIF
437C--------------------
438C UPDATE YIELD
439C--------------------
440 300 CONTINUE
441 DO i=1,nel
442 IF(alpha(i) > one) THEN
443 IF(ipla(i) == 1) THEN
444 IF(pnew(i) < ps) THEN
445 alphayld(i)= one+ eta*(ps-pnew(i))**en
446 ELSE
447 alphayld(i)=one
448 ENDIF
449 ENDIF
450 ENDIF
451 ENDDO
452C--------------------
453C DEVIATORIC STRESSES
454C--------------------
455 IF(iflag2==2)THEN
456 DO i=1,nel
457 dav = (depsxx(i)+depsyy(i)+depszz(i))*third
458 g = max(zero,gs+dgdalpha*(alpha(i)-one))
459 g2 = two*g
460 signxx(i)=sigoxx(i)+pold(i)+g2*(depsxx(i)-dav)
461 signyy(i)=sigoyy(i)+pold(i)+g2*(depsyy(i)-dav)
462 signzz(i)=sigozz(i)+pold(i)+g2*(depszz(i)-dav)
463 signxy(i)=sigoxy(i)+g*depsxy(i)
464 signyz(i)=sigoyz(i)+g*depsyz(i)
465 signzx(i)=sigozx(i)+g*depszx(i)
466 ENDDO
467 ENDIF
468C-------
469C UPDATE
470C-------
471 DO i=1,nel
472 uvar(i,2)=theta(i)
473 uvar(i,3)=alpha(i)
474 uvar(i,4)=alphayld(i)
475 signxx(i)=signxx(i)-pnew(i)
476 signyy(i)=signyy(i)-pnew(i)
477 signzz(i)=signzz(i)-pnew(i)
478C SET MAXIMUM VISCOSITY
479 viscmax(i)=zero
480 ENDDO
481C SET SOUND SPEED
482 DO i=1,nel
483 IF(alpha(i) > one) THEN
484 soundsp(i)=sqrt(bulk0/rhos0(1))*(aa*(alpha(i)-one)+one)
485 ELSE
486 soundsp(i)=sqrt(bulk0/rhos0(1))
487 ENDIF
488 ENDDO
489C
490 RETURN
my_real
#define my_real
Definition cppsort.cpp:32
alpha
#define alpha
Definition eval.h:35
max
#define max(a, b)
Definition macros.h:21
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
message_mod
Definition message_mod.F:1249
ancmsg
subroutine ancmsg(msgid, msgtype, anmode, i1, i2, i3, i4, i5, i6, i7, i8, i9, i10, i11, i12, i13, i14, i15, i16, i17, i18, i19, i20, r1, r2, r3, r4, r5, r6, r7, r8, r9, c1, c2, c3, c4, c5, c6, c7, c8, c9, prmode)
Definition message.F:889
arret
subroutine arret(nn)
Definition arret.F:87