OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps88.F File Reference
#include "implicit_f.inc"
#include "scr05_c.inc"
#include "impl1_c.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine sigeps88 (nel, nuparam, nuvar, nfunc, ifunc, npf, tf, time, timestep, uparam, rho0, rho, volume, eint, ngl, 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, ismstr, israte, asrate, et, ihet, offg, epsth3, iexpan, epsd)

Function/Subroutine Documentation

◆ sigeps88()

subroutine sigeps88 ( integer nel,
integer nuparam,
integer nuvar,
integer nfunc,
integer, dimension(nfunc) ifunc,
integer, dimension(*) npf,
tf,
time,
timestep,
uparam,
rho0,
rho,
volume,
eint,
integer, dimension(*) ngl,
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,
integer ismstr,
integer israte,
asrate,
et,
integer ihet,
offg,
epsth3,
integer iexpan,
epsd )

Definition at line 32 of file sigeps88.F.

45C-----------------------------------------------
46C I M P L I C I T T Y P E S
47C-----------------------------------------------
48#include "implicit_f.inc"
49C-----------------------------------------------
50C C O M M O N
51C-----------------------------------------------
52#include "scr05_c.inc"
53#include "impl1_c.inc"
54C----------------------------------------------------------------
55C I N P U T A R G U M E N T S
56C----------------------------------------------------------------
57 INTEGER NEL, NUPARAM, NUVAR,ISMSTR,NGL(*),IHET,
58 . ISRATE,IEXPAN
59 my_real
60 . time , timestep , uparam(nuparam),
61 . rho(nel), rho0(nel), volume(nel), eint(nel),
62 . epspxx(nel), epspyy(nel), epspzz(nel),
63 . epspxy(nel), epspyz(nel), epspzx(nel),
64 . depsxx(nel), depsyy(nel), depszz(nel),
65 . depsxy(nel), depsyz(nel), depszx(nel),
66 . epsxx(nel), epsyy(nel), epszz(nel),
67 . epsxy(nel), epsyz(nel), epszx(nel),
68 . sigoxx(nel), sigoyy(nel), sigozz(nel),
69 . sigoxy(nel), sigoyz(nel), sigozx(nel),
70 . asrate,offg(nel),epsth3(nel),epsd(nel)
71C----------------------------------------------------------------
72C O U T P U T A R G U M E N T S
73C----------------------------------------------------------------
74 my_real
75 . signxx(nel), signyy(nel), signzz(nel),
76 . signxy(nel), signyz(nel), signzx(nel),
77 . sigvxx(nel), sigvyy(nel), sigvzz(nel),
78 . sigvxy(nel), sigvyz(nel), sigvzx(nel),
79 . soundsp(nel), viscmax(nel), et(nel)
80C----------------------------------------------------------------
81C I N P U T O U T P U T A R G U M E N T S
82C----------------------------------------------------------------
83 my_real
84 . uvar(nel,nuvar), off(nel)
85C----------------------------------------------------------------
86C VARIABLES FOR FUNCTION INTERPOLATION
87C----------------------------------------------------------------
88 INTEGER NPF(*), NFUNC, IFUNC(NFUNC)
89 my_real finter,fintte,tf(*),fint2v
90 EXTERNAL finter,fintte
91C----------------------------------------------------------------
92C L O C A L V A R I B L E S
93C----------------------------------------------------------------
94 INTEGER I,J,ITENS,NLOAD, NUNLOAD,NE_LOAD,NE_UNLOAD,IUNLOAD,
95 . JJ,J1,J2,IUNL_FOR,KK,NN,K,ICASE,IDIR,
96 . INDX_L(NEL),INDX_UNL(NEL),ITYPE,IFD(3),ITAG(NEL)
97 my_real
98 . av(6,nel),ev(3,nel),ee(3,nel),edot(3,nel),
99 . evv(3,nel),dirprv(3,3,nel),erate(nel),rv(nel),
100 . yfac(nfunc),rate(nfunc),t1(nel),t2(nel),t3(nel),
101 . ee0(6,nel),ee1(6,nel),fg1(6),gmax(nel),fqstat(nel,3),
102 . dam(nel),funl(3),
103 . deint,p,invr,df,fac,d,yfac_unl,f01,f02,f03,eps1,
104 . rbulk,nu,gs,hys,shape,f(3),f1,f2,dtinv,
105 . pvisc,dd,beta,dd1,dd2,dd3,aa,
106 . de1,de2,de3,de4,de5,de6, ert11,ert12,ert13,ert21,ert22,ert23,
107 . ert31,ert32,ert33,dep1,dep2,dep3,df1,df2,df3, xint,yint, lam_tens,
108 . t_tens, scale(3),es(3),evs(3),fstar,e2,f3,f_int,xinc,yinc,
109 . lam_comp,t_comp,e_tens,e_comp,yintl,bb,fc,ft,ec_start,et_start,
110 . ecurent(nel),epsp(3,nel),epe,dmax,sum,ft1,ft2,fc1,fc2,fcd,ftd,fcs,fts,
111 . epsdold(nel),yint0,yinc0,fts_unl,fcs_unl,dmoy,eintv,epseq(nel),
112 . t_tens_unl, t_comp_unl,scale_unl,scale_l,dam_f(nel),lam_xinc,
113 . lam_xint
114C----------------------------------------------------------------
115C material parameters
116 rbulk = uparam(1)
117 nu = uparam(2)
118 gs = uparam(3)
119 nload = int(uparam(4))
120 iunl_for = int(uparam(5))
121 hys = uparam(6)
122 shape = uparam(7)
123 icase = nint(uparam(9))
124C
125 beta = nu
126C
127 DO j=1,nload
128 rate(j) = uparam(9 + 2*j - 1 )
129 yfac(j) = uparam(9 + 2*j )
130 ENDDO
131 yfac_unl = uparam(9 + 2*nload + 2 )
132C iniialisation des variabesl users
133 IF(time == zero)THEN
134 DO j = 1,nuvar
135 DO i = 1, nel
136 uvar(i,j) = zero
137 ENDDO
138 ENDDO
139 DO i = 1, nel
140 uvar(i,22) = one
141 uvar(i,25) = one
142 uvar(i,30) = zero ! -0.99
143 uvar(i,31) = zero ! ONE
144 uvar(i,18) = one
145 uvar(i,21) = zero
146 uvar(i,24) = zero
147 ENDDO
148 ENDIF
149C
150 DO i=1,nel
151 av(1,i)=epsxx(i)
152 av(2,i)=epsyy(i)
153 av(3,i)=epszz(i)
154 av(4,i)=half*epsxy(i)
155 av(5,i)=half*epsyz(i)
156 av(6,i)=half*epszx(i)
157 ENDDO
158
159CEigenvalues needed to be calculated in double precision
160C for a simple precision executing
161 IF (iresp==1) THEN
162 CALL valpvecdp(av,evv,dirprv,nel)
163 ELSE
164 CALL valpvec(av,evv,dirprv,nel)
165 ENDIF
166C-ISMSTR=0-NO SMALL STRAIN OPTION:STRAINS ARE LOGARITHMIC, STRESS IS CAUCHY
167C-ISMSTR=1-SMALL STRAIN OPTION:STRAINS ARE ENGINEERING, STRESS IS CAUCHY
168C-ISMSTR=2-SMALL STRAIN OPTION:STRAINS ARE ENGINEERING, STRESS IS BIOT
169C-ISMSTR=3-NO SMALL STRAIN OPTION:STRESS IS BIOT
170 IF(ismstr==0.OR.ismstr==2.OR.ismstr==4) THEN
171 DO i=1,nel
172C ---- (STRAIN IS LOGARITHMIC)
173 ev(1,i)=exp(evv(1,i))
174 ev(2,i)=exp(evv(2,i))
175 ev(3,i)=exp(evv(3,i))
176 ENDDO
177 ELSEIF(ismstr==10.OR.ismstr==12) THEN
178 DO i =1,nel
179 IF(offg(i)<=one) THEN
180 ev(1,i)=sqrt(evv(1,i)+ one)
181 ev(2,i)=sqrt(evv(2,i)+ one)
182 ev(3,i)=sqrt(evv(3,i)+ one)
183 ELSE
184 ev(1,i)=evv(1,i)+ one
185 ev(2,i)=evv(2,i)+ one
186 ev(3,i)=evv(3,i)+ one
187 END IF
188 ENDDO
189 ELSE
190C ---- STRAIN IS ENGINEERING)
191 DO i=1,nel
192 ev(1,i)=evv(1,i) + one
193 ev(2,i)=evv(2,i) + one
194 ev(3,i)=evv(3,i) + one
195 ENDDO
196 ENDIF
197
198 DO i=1,nel
199 ee(1,i) = ev(1,i) - one
200 ee(2,i) = ev(2,i) - one
201 ee(3,i) = ev(3,i) - one
202C from principal to global
203 ee1(1,i) = dirprv(1,1,i)*dirprv(1,1,i)*ee(1,i)
204 . + dirprv(1,2,i)*dirprv(1,2,i)*ee(2,i)
205 . + dirprv(1,3,i)*dirprv(1,3,i)*ee(3,i)
206
207 ee1(2,i) = dirprv(2,2,i)*dirprv(2,2,i)*ee(2,i)
208 . + dirprv(2,3,i)*dirprv(2,3,i)*ee(3,i)
209 . + dirprv(2,1,i)*dirprv(2,1,i)*ee(1,i)
210
211 ee1(3,i) = dirprv(3,3,i)*dirprv(3,3,i)*ee(3,i)
212 . + dirprv(3,1,i)*dirprv(3,1,i)*ee(1,i)
213 . + dirprv(3,2,i)*dirprv(3,2,i)*ee(2,i)
214
215 ee1(4,i) = dirprv(1,1,i)*dirprv(2,1,i)*ee(1,i)
216 . + dirprv(1,2,i)*dirprv(2,2,i)*ee(2,i)
217 . + dirprv(1,3,i)*dirprv(2,3,i)*ee(3,i)
218
219 ee1(5,i) = dirprv(2,2,i)*dirprv(3,2,i)*ee(2,i)
220 . + dirprv(2,3,i)*dirprv(3,3,i)*ee(3,i)
221 . + dirprv(2,1,i)*dirprv(3,1,i)*ee(1,i)
222
223 ee1(6,i) = dirprv(3,3,i)*dirprv(1,3,i)*ee(3,i)
224 . + dirprv(3,1,i)*dirprv(1,1,i)*ee(1,i)
225 . + dirprv(3,2,i)*dirprv(1,2,i)*ee(2,i)
226
227 ENDDO
228C Strain rate by principal direction.
229 dtinv =timestep/max(timestep**2,em20)
230C
231 DO i= 1,nel
232C
233 ee0(1,i) = uvar(i,1)
234 ee0(2,i) = uvar(i,2)
235 ee0(3,i) = uvar(i,3)
236 ee0(4,i) = uvar(i,4)
237 ee0(5,i) = uvar(i,5)
238 ee0(6,i) = uvar(i,6)
239C
240 de1 = ee1(1,i) - ee0(1,i)
241 de2 = ee1(2,i) - ee0(2,i)
242 de3 = ee1(3,i) - ee0(3,i)
243 de4 = ee1(4,i) - ee0(4,i)
244 de5 = ee1(5,i) - ee0(5,i)
245 de6 = ee1(6,i) - ee0(6,i)
246
247 ert11 =dirprv(1,1,i)*de1 + dirprv(2,1,i)*de4 + dirprv(3,1,i)*de6
248 ert12 =dirprv(1,2,i)*de1 + dirprv(2,2,i)*de4 + dirprv(3,2,i)*de6
249 ert13 =dirprv(1,3,i)*de1 + dirprv(2,3,i)*de4 + dirprv(3,3,i)*de6
250
251 ert21 =dirprv(1,1,i)*de4 + dirprv(2,1,i)*de2 + dirprv(3,1,i)*de5
252 ert22 =dirprv(1,2,i)*de4 + dirprv(2,2,i)*de2 + dirprv(3,2,i)*de5
253 ert23 =dirprv(1,3,i)*de4 + dirprv(2,3,i)*de2 + dirprv(3,3,i)*de5
254
255 ert31 =dirprv(1,1,i)*de6 + dirprv(2,1,i)*de5 + dirprv(3,1,i)*de3
256 ert32 =dirprv(1,2,i)*de6 + dirprv(2,2,i)*de5 + dirprv(3,2,i)*de3
257 ert33 =dirprv(1,3,i)*de6 + dirprv(2,3,i)*de5 + dirprv(3,3,i)*de3
258C
259 dep1 = dirprv(1,1,i)*ert11 + dirprv(2,1,i)*ert21
260 . + dirprv(3,1,i)*ert31
261 dep2 = dirprv(1,2,i)*ert12 + dirprv(2,2,i)*ert22
262 . + dirprv(3,2,i)*ert32
263 dep3 = dirprv(1,3,i)*ert13 + dirprv(2,3,i)*ert23
264 . + dirprv(3,3,i)*ert33
265C
266 edot(1,i) = abs(dep1)*dtinv
267 edot(2,i) = abs(dep2)*dtinv
268 edot(3,i) = abs(dep3)*dtinv
269C
270 erate(i) = max(edot(1,i),edot(2,i),edot(3,i))
271 rv(i) = ev(1,i)*ev(2,i)*ev(3,i)
272 ENDDO
273C----THERM STRESS COMPUTATION-----
274 IF(iexpan > 0.AND.(ismstr==10.OR.ismstr==11.OR.ismstr==12)) THEN
275 DO i=1,nel
276 rv(i) = rv(i) - epsth3(i)
277 ENDDO
278 ENDIF
279C----------------
280 IF(israte > 0) THEN
281 DO i=1,nel
282 edot(1,i) = asrate*edot(1,i) + (one - asrate)*uvar(i,7)
283 edot(2,i) = asrate*edot(2,i) + (one - asrate)*uvar(i,8)
284 edot(3,i) = asrate*edot(3,i) + (one - asrate)*uvar(i,9)
285 uvar(i,7) = edot(1,i)
286 uvar(i,8) = edot(2,i)
287 uvar(i,9) = edot(3,i)
288 epsdold(i) = epsd(i)
289 epsd(i) = max(edot(1,i),edot(2,i),edot(3,i))
290 edot(1,i) = epsd(i)
291 edot(2,i) = epsd(i)
292 edot(3,i) = epsd(i)
293 ENDDO
294 ELSE
295C for outp only
296 DO i=1,nel
297 uvar(i,7) = edot(1,i)
298 uvar(i,8) = edot(2,i)
299 uvar(i,9) = edot(3,i)
300 epsdold(i) = epsd(i)
301 epsd(i) = max(edot(1,i),edot(2,i),edot(3,i))
302 edot(1,i) = epsd(i)
303 edot(2,i) = epsd(i)
304 edot(3,i) = epsd(i)
305 ENDDO
306 ENDIF
307C----------------
308C----------------
309C compute total energy
310C ---------------
311 ne_load = 0
312 ne_unload = 0
313 DO i=1,nel
314C from global frame to principal
315C
316 f(1) = ev(1,i)*yfac(1)*finter(ifunc(1),ee(1,i),npf,tf,df)
317 f(2) = ev(2,i)*yfac(1)*finter(ifunc(1),ee(2,i),npf,tf,df)
318 f(3) = ev(3,i)*yfac(1)*finter(ifunc(1),ee(3,i),npf,tf,df)
319 gmax(i) = gs
320 epseq(i) = sqrt(ee(1,i)**2 + ee(2,i)**2 + ee(3,i)**2)
321 fac = one
322 invr = one/max(em20,rv(i))
323 p = invr*rbulk*(rv(i) - one)
324 t1(i) =invr*(two_third*f(1) - third*(f(2) + f(3))) + p
325 t2(i) =invr*(two_third*f(2) - third*(f(1) + f(3))) + p
326 t3(i) =invr*(two_third*f(3) - third*(f(1) + f(2))) + p
327 ecurent(i)= half*(ee(1,i)*t1(i) + ee(2,i)*t2(i) + ee(3,i)*t3(i))
328 ecurent(i) = max(em20,ecurent(i))
329 IF(ecurent(i) >= uvar(i,17)) uvar(i,20) = zero
330 uvar(i,17)= max(uvar(i,17), ecurent(i))
331 deint = ecurent(i) - uvar(i,16)
332 uvar(i,16) = ecurent(i)
333 f(1) = f(1)/ev(1,i)
334 f(2) = f(2)/ev(2,i)
335 f(3) = f(3)/ev(3,i)
336C
337 IF(uvar(i,19) == -one) THEN
338 IF(deint/uvar(i,16) >= 1e-7) THEN
339 ne_load = ne_load + 1
340 indx_l(ne_load) = i
341 uvar(i,19) = one
342 uvar(i,20) = one ! re-loading
343C ....
344 IF(iunl_for > 1) THEN
345 uvar(i,22) = one
346 uvar(i,25) = one
347 uvar(i,26) = zero
348 uvar(i,27) = one
349 uvar(i,28) = zero
350 uvar(i,29) = one
351 ENDIF !
352 ELSE
353 ne_unload = ne_unload + 1
354 indx_unl(ne_unload) = i
355 uvar(i,19) = -one
356 ENDIF
357 ELSE
358 IF(deint/max(em20,uvar(i,17)) >= zero) THEN
359 ne_load = ne_load + 1
360 indx_l(ne_load) = i
361 uvar(i,19) = one
362 ELSE
363 ne_unload = ne_unload + 1
364 indx_unl(ne_unload) = i
365 uvar(i,19) = - one
366! the curve must be monotonic
367 IF(nload > 1 ) THEN
368 lam_tens = max(one,ev(1,i),ev(2,i),ev(3,i))
369 e_tens = max(zero,ee(1,i),ee(2,i),ee(3,i))
370 t_tens= max(zero,f(1),f(2),f(3))
371 lam_comp = min(one,ev(1,i),ev(2,i),ev(3,i))
372 e_comp = min(zero,ee(1,i),ee(2,i),ee(3,i))
373 t_comp = min(zero,f(1),f(2),f(3))
374 IF(iunl_for == 1) THEN
375 uvar(i,22) = one
376 uvar(i,25) = one
377 IF(uvar(i,20) == zero ) THEN
378 kk = nload + 1
379 j1 = 1
380 DO k=2,nload - 1
381 IF( epsdold(i) >= rate(k) )j1 = k
382 ENDDO
383 j2 = j1 + 1
384 aa = (epsdold(i) - rate(j1))/(rate(j2) - rate(j1))
385 df = zero
386C
387 ft1 =lam_tens*yfac(j1)*finter(ifunc(j1),e_tens,npf,tf,df1)
388 ft2 =lam_tens*yfac(j2)*finter(ifunc(j2),e_tens,npf,tf,df2)
389 fts =lam_tens*yfac_unl*finter(ifunc(kk),e_tens,npf,tf,df)
390C compression
391 fc1 =lam_comp*yfac(j1)*finter(ifunc(j1),e_comp,npf,tf,df1)
392 fc2 =lam_comp*yfac(j2)*finter(ifunc(j2),e_comp,npf,tf,df2)
393 fcs =lam_comp*yfac_unl*finter(ifunc(kk),e_comp,npf,tf,df1)
394 ftd = ft1 + aa*(ft2 - ft1)
395 fcd = fc1 + aa*(fc2 - fc1)
396 uvar(i,22) = ftd/fts
397 uvar(i,25) = fcd/fcs
398 uvar(i,22)= uvar(i,22)/e_tens
399 uvar(i,25)= uvar(i,25)/abs(e_comp)
400 ELSE ! iunl_for > 1
401 j1 = 1
402 DO k=2,nload - 1
403 IF( epsdold(i) >= rate(k) )j1 = k
404 ENDDO
405 j2 = j1 + 1
406 aa = (epsdold(i) - rate(j1))/(rate(j2) - rate(j1))
407 df = zero
408C
409 ft1 =lam_tens*yfac(j1)*finter(ifunc(j1),e_tens,npf,tf,df1)
410 ft2 =lam_tens*yfac(j2)*finter(ifunc(j2),e_tens,npf,tf,df2)
411 fts =lam_tens*yfac(1)*finter(ifunc(1),e_tens,npf,tf,df)
412C compression
413 fc1 =lam_comp*yfac(j1)*finter(ifunc(j1),e_comp,npf,tf,df1)
414 fc2 =lam_comp*yfac(j2)*finter(ifunc(j2),e_comp,npf,tf,df2)
415 fcs =lam_comp*yfac(1)*finter(ifunc(1),e_comp,npf,tf,df1)
416 fac = one
417 DO nn=1,20
418 fac = fac*(-half)
419 dd1 = lam_tens**fac - one
420 ft1 =ft1 + (dd1 + one)*yfac(j1)*finter(ifunc(j1),dd1,npf,tf,df)
421 ft2 =ft2 + (dd1 + one)*yfac(j2)*finter(ifunc(j2),dd1,npf,tf,df)
422 fts =fts + (dd1 + one)*yfac(1)*finter(ifunc(1),dd1,npf,tf,df)
423 dd1 = lam_comp**fac - one
424 fc1 =fc1 + (dd1 + one)*yfac(j1)*finter(ifunc(j1),dd1,npf,tf,df)
425 fc2 =fc2 + (dd1 + one)*yfac(j2)*finter(ifunc(j2),dd1,npf,tf,df)
426 fcs =fcs + (dd1 + one)*yfac(1)*finter(ifunc(1),dd1,npf,tf,df)
427 ENDDO ! NN
428 ftd = ft1 + aa*(ft2 - ft1)
429 fcd = fc1 + aa*(fc2 - fc1)
430 uvar(i,22) = (ftd/fts - one)/e_tens ! tension
431 uvar(i,25) = (fcd/fcs - one)/e_comp ! compression
432 ENDIF ! UVAR 20
433 ENDIF ! i IFOR_UNL
434 ENDIF ! strain rate effet !NLOAD > 1
435 ENDIF ! DEINT
436 ENDIF ! UVAR 19
437 uvar(i,32)= epseq(i)
438 ENDDO
439C
440 IF(ne_load > 0 ) THEN
441 IF(nload == 1) THEN
442 IF(iunl_for == 1) THEN
443 kk = nload + 1
444 DO jj=1,ne_load
445 i = indx_l(jj)
446 invr = one/max(em20,rv(i))
447 dam(i) = uvar(i,16)/uvar(i,17)
448 !!
449 f(1) = ev(1,i)*yfac(1)*finter(ifunc(1),ee(1,i),npf,tf,df1)
450 f(2) = ev(2,i)*yfac(1)*finter(ifunc(1),ee(2,i),npf,tf,df2)
451 f(3) = ev(3,i)*yfac(1)*finter(ifunc(1),ee(3,i),npf,tf,df3)
452 gmax(i) = max(gmax(i),yfac(1)*df1,yfac(1)*df2,
453 . yfac(1)*df3 )
454 IF(dam(i) /= one) THEN
455 funl(1)=ev(1,i)*yfac_unl*finter(ifunc(kk),ee(1,i),npf,tf,df)
456 funl(2)=ev(2,i)*yfac_unl*finter(ifunc(kk),ee(2,i),npf,tf,df)
457 funl(3)=ev(3,i)*yfac_unl*finter(ifunc(kk),ee(3,i),npf,tf,df)
458 !!
459 IF(f(1) /= zero) dam_f(i) = max(dam_f(i),funl(1)/f(1))
460 IF(f(2) /= zero) dam_f(i) = max(dam_f(i),funl(2)/f(2))
461 IF(f(3) /= zero) dam_f(i) = max(dam_f(i),funl(3)/f(3))
462 IF(dam_f(i) <= one)THEN
463 dam(i) = max(dam(i),dam_f(i))
464 ELSE
465 dam(i) = one
466 ENDIF
467 ENDIF
468 dam(i) = max(dam(i),uvar(i,18))
469 dam(i) = min(one, dam(i))
470 fac = one
471 DO nn=1,20
472 fac = fac*(-half)
473 dd1 = ev(1,i)**fac - one
474 dd2 = ev(2,i)**fac - one
475 dd3 = ev(3,i)**fac - one
476 f(1) = f(1) +
477 . (dd1 + one)*yfac(1)*finter(ifunc(1),dd1,npf,tf,df)
478 f(2) = f(2) +
479 . (dd2 + one)*yfac(1)*finter(ifunc(1),dd2,npf,tf,df)
480 f(3) = f(3) +
481 . (dd3 + one)*yfac(1)*finter(ifunc(1),dd3,npf,tf,df)
482 ENDDO
483 !
484 p = invr*rbulk*(rv(i) - one)
485 t1(i) = dam(i)*(invr*(two_third*f(1) - third*(f(2) + f(3))) )+ p
486 t2(i) = dam(i)*(invr*(two_third*f(2) - third*(f(1) + f(3))) )+ p
487 t3(i) = dam(i)*(invr*(two_third*f(3) - third*(f(1) + f(2))) )+ p
488 uvar(i,19) = one
489 ENDDO ! NE_LOAD
490 ELSE ! NLOAD > 0 ! with strain rate effect
491 DO jj=1,ne_load
492 i = indx_l(jj)
493 invr = one/max(em20,rv(i))
494 f(1) = ev(1,i)*yfac(1)*finter(ifunc(1),ee(1,i),npf,tf,df1)
495 f(2) = ev(2,i)*yfac(1)*finter(ifunc(1),ee(2,i),npf,tf,df2)
496 f(3) = ev(3,i)*yfac(1)*finter(ifunc(1),ee(3,i),npf,tf,df3)
497 gmax(i) = max(gmax(i),yfac(1)*df1,yfac(1)*df2,
498 . yfac(1)*df3 )
499 fac = one
500 DO nn=1,20
501 fac = fac*(-half)
502 dd1 = ev(1,i)**fac - one
503 dd2 = ev(2,i)**fac - one
504 dd3 = ev(3,i)**fac - one
505 f(1) = f(1) +
506 . (dd1 + one)*yfac(1)*finter(ifunc(1),dd1,npf,tf,df)
507 f(2) = f(2) +
508 . (dd2 + one)*yfac(1)*finter(ifunc(1),dd2,npf,tf,df)
509 f(3) = f(3) +
510 . (dd3 + one)*yfac(1)*finter(ifunc(1),dd3,npf,tf,df)
511 ENDDO
512 p = invr*rbulk*(rv(i) - one)
513 t1(i) =uvar(i,22)*invr*(two_third*f(1) - third*(f(2) + f(3))) + p
514 t2(i) =uvar(i,22)*invr*(two_third*f(2) - third*(f(1) + f(3))) + p
515 t3(i) =uvar(i,22)*invr*(two_third*f(3) - third*(f(1) + f(2))) + p
516 uvar(i,19) = one
517 ENDDO
518 ENDIF ! iunl_for
519 ELSE ! NLOAD > 1
520 IF(iunl_for == 1) THEN
521 kk = nload + 1
522 DO jj=1,ne_load
523 i = indx_l(jj)
524 IF(uvar(i,20) == one .AND. ecurent(i) <= uvar(i,17))THEN
525 f(1)=ev(1,i)*yfac_unl*finter(ifunc(kk),ee(1,i),npf,tf,df)
526 f(2)=ev(2,i)*yfac_unl*finter(ifunc(kk),ee(2,i),npf,tf,df)
527 f(3)=ev(3,i)*yfac_unl*finter(ifunc(kk),ee(3,i),npf,tf,df)
528 e_tens = max(zero,ee(1,i),ee(2,i),ee(3,i))
529 dam(i)= uvar(i,22)*e_tens
530 IF(itype == -1 .OR. (itype == 2 .AND. rv(i) < one)) THEN
531 e_comp = min(zero,ee(1,i),ee(2,i),ee(3,i))
532 dam(i)= uvar(i,25)*abs(e_comp)
533 ENDIF
534 fac = one
535 DO nn=1,20
536 fac = (-half)*fac
537 dd1 = ev(1,i)**fac - one
538 dd2 = ev(2,i)**fac - one
539 dd3 = ev(3,i)**fac - one
540 f(1) = f(1) +
541 . (dd1 + one)*yfac_unl*finter(ifunc(kk),dd1,npf,tf,df)
542 f(2) = f(2) +
543 . (dd2 + one)*yfac_unl*finter(ifunc(kk),dd2,npf,tf,df)
544 f(3) = f(3) +
545 . (dd3 + one)*yfac_unl*finter(ifunc(kk),dd3,npf,tf,df)
546 ENDDO
547 invr = one/max(em20,rv(i))
548 p = invr*rbulk*(rv(i) - one)
549 t1(i) = invr*dam(i)*(two_third*f(1) - third*(f(2) + f(3))) + p
550 t2(i) = invr*dam(i)*(two_third*f(2) - third*(f(1) + f(3))) + p
551 t3(i) = invr*dam(i)*(two_third*f(3) - third*(f(1) + f(2))) + p
552 ELSE
553 DO j=1,3
554 j1 = 1
555 DO k=2,nload - 1
556 IF( edot(j,i) >= rate(k) )j1 = k
557 ENDDO
558 j2 = j1 + 1
559 aa = (edot(j,i) - rate(j1))/(rate(j2) - rate(j1))
560 df = zero
561C
562 f1 =ev(j,i)*yfac(j1)*finter(ifunc(j1),ee(j,i),npf,tf,df1)
563 f2 =ev(j,i)*yfac(j2)*finter(ifunc(j2),ee(j,i),npf,tf,df2)
564 df1 = ev(j,i)*df1
565 df2 = ev(j,i)*df2
566 df = max(df, df1 + aa*(df2 - df1))
567 gmax(i) = max(gmax(i),yfac(j1)*df1,yfac(j2)*df2 )
568 fac = one
569 DO nn=1,20
570 fac = fac*(-half)
571 dd1 = ev(j,i)**fac - one
572 f1 =f1+(dd1 + one)*yfac(j1)*finter(ifunc(j1),dd1,npf,tf,df)
573 f2 =f2+(dd1 + one)*yfac(j2)*finter(ifunc(j2),dd1,npf,tf,df)
574 ENDDO ! NN
575 f(j) = f1 + aa*(f2 - f1)
576 ENDDO ! J
577 invr = one/max(em20,rv(i))
578 p = invr*rbulk*(rv(i) - one)
579 t1(i) =invr*(two_third*f(1) - third*(f(2) + f(3))) + p
580 t2(i) =invr*(two_third*f(2) - third*(f(1) + f(3))) + p
581 t3(i) =invr*(two_third*f(3) - third*(f(1) + f(2))) + p
582 ENDIF
583 ENDDO
584 ElSE ! iunl_for > 1
585 DO jj=1,ne_load
586 i = indx_l(jj)
587 DO j=1,3
588 j1 = 1
589 DO k=2,nload - 1
590 IF( edot(j,i) >= rate(k) )j1 = k
591 ENDDO
592 j2 = j1 + 1
593 aa = (edot(j,i) - rate(j1))/(rate(j2) - rate(j1))
594 df = zero
595C
596 f1 =ev(j,i)*yfac(j1)*finter(ifunc(j1),ee(j,i),npf,tf,df1)
597 f2 =ev(j,i)*yfac(j2)*finter(ifunc(j2),ee(j,i),npf,tf,df2)
598 df1 = ev(j,i)*df1
599 df2 = ev(j,i)*df2
600 df = max(df, df1 + aa*(df2 - df1))
601 gmax(i) = max(gmax(i),yfac(j1)*df1,yfac(j2)*df2 )
602 fac = one
603 DO nn=1,20
604 fac = fac*(-half)
605 dd1 = ev(j,i)**fac - one
606 f1 =f1+(dd1 + one)*yfac(j1)*finter(ifunc(j1),dd1,npf,tf,df)
607 f2 =f2+(dd1 + one)*yfac(j2)*finter(ifunc(j2),dd1,npf,tf,df)
608 ENDDO ! NN
609 f(j) = f1 + aa*(f2 - f1)
610 ENDDO ! J
611 invr = one/max(em20,rv(i))
612 p = invr*rbulk*(rv(i) - one)
613 t1(i) =invr*(two_third*f(1) - third*(f(2) + f(3))) + p
614 t2(i) =invr*(two_third*f(2) - third*(f(1) + f(3))) + p
615 t3(i) =invr*(two_third*f(3) - third*(f(1) + f(2))) + p
616 ENDDO
617 ENDIF ! IUNL_FOR
618 ENDIF ! NLOAD
619 ENDIF ! NE_LOAD
620
621C unloading
622C
623 IF(ne_unload > 0 ) THEN
624C
625 SELECT CASE (iunl_for)
626 !
627 CASE (0) ! Nonlinear elastic behavior (with out hysteresis)
628 DO jj=1,ne_unload
629 i = indx_unl(jj)
630 f(1) = ev(1,i)*yfac(1)*finter(ifunc(1),ee(1,i),npf,tf,df1)
631 f(2) = ev(2,i)*yfac(1)*finter(ifunc(1),ee(2,i),npf,tf,df2)
632 f(3) = ev(3,i)*yfac(1)*finter(ifunc(1),ee(3,i),npf,tf,df3)
633 gmax(i) = max(gmax(i),yfac(1)*df1,yfac(1)*df2,
634 . yfac(1)*df3 )
635 fac = one
636 DO nn=1,20
637 fac = fac*(-half)
638 dd1 = ev(1,i)**fac - one
639 dd2 = ev(2,i)**fac - one
640 dd3 = ev(3,i)**fac - one
641 f(1) = f(1) +
642 . (dd1 + one)*yfac(1)*finter(ifunc(1),dd1,npf,tf,df)
643 f(2) = f(2) +
644 . (dd2 + one)*yfac(1)*finter(ifunc(1),dd2,npf,tf,df)
645 f(3) = f(3) +
646 . (dd3 + one)*yfac(1)*finter(ifunc(1),dd3,npf,tf,df)
647 ENDDO
648 invr = one/max(em20,rv(i))
649 p = invr*rbulk*(rv(i) - one)
650 t1(i) = invr*(two_third*f(1) - third*(f(2) + f(3))) + p
651 t2(i) = invr*(two_third*f(2) - third*(f(1) + f(3))) + p
652 t3(i) = invr*(two_third*f(3) - third*(f(1) + f(2))) + p
653 ENDDO ! NE_UNLOAD
654C
655 CASE (1) ! unloading following curve
656C
657 IF(nload == 1) THEN
658 kk = nload + 1
659 DO jj=1,ne_unload
660 i = indx_unl(jj)
661C
662 dam(i) = uvar(i,16)/uvar(i,17)
663 f(1) = ev(1,i)*yfac(1)*finter(ifunc(1),ee(1,i),npf,tf,df1)
664 f(2) = ev(2,i)*yfac(1)*finter(ifunc(1),ee(2,i),npf,tf,df2)
665 f(3) = ev(3,i)*yfac(1)*finter(ifunc(1),ee(3,i),npf,tf,df3)
666 gmax(i) = max(gmax(i),yfac(1)*df1,yfac(1)*df2,
667 . yfac(1)*df3 )
668 funl(1)=ev(1,i)*yfac_unl*finter(ifunc(kk),ee(1,i),npf,tf,df)
669 funl(2)=ev(2,i)*yfac_unl*finter(ifunc(kk),ee(2,i),npf,tf,df)
670 funl(3)=ev(3,i)*yfac_unl*finter(ifunc(kk),ee(3,i),npf,tf,df)
671C
672 IF(f(1)/=zero) dam_f(i) = max(dam_f(i),funl(1)/f(1))
673 IF(f(2)/=zero) dam_f(i) = max(dam_f(i),funl(2)/f(2))
674 IF(f(3)/=zero) dam_f(i) = max(dam_f(i),funl(3)/f(3))
675 IF(dam_f(i) <= one)THEN
676 dam(i) = max(dam(i),dam_f(i))
677 ELSE
678 dam(i) = one
679 ENDIF
680 dam(i) = min(one,dam(i))
681 uvar(i,18) = dam(i)
682 fac = one
683 DO nn=1,20
684 fac = fac*(-half)
685 dd1 = ev(1,i)**fac - one
686 dd2 = ev(2,i)**fac - one
687 dd3 = ev(3,i)**fac - one
688 f(1) = f(1) +
689 . (dd1 + one)*yfac(1)*finter(ifunc(1),dd1,npf,tf,df)
690 f(2) = f(2) +
691 . (dd2 + one)*yfac(1)*finter(ifunc(1),dd2,npf,tf,df)
692 f(3) = f(3) +
693 . (dd3 + one)*yfac(1)*finter(ifunc(1),dd3,npf,tf,df)
694 ENDDO
695 invr = one/max(em20,rv(i))
696 p = invr*rbulk*(rv(i) - one)
697 t1(i) = dam(i)*(invr*(two_third*f(1) - third*(f(2) + f(3))) )+ p
698 t2(i) = dam(i)*(invr*(two_third*f(2) - third*(f(1) + f(3))) )+ p
699 t3(i) = dam(i)*(invr*(two_third*f(3) - third*(f(1) + f(2))) )+ p
700 ENDDO ! NE_UNLOAD
701 ELSE
702 kk = nload + 1
703 DO jj=1,ne_unload
704 i = indx_unl(jj)
705 dam(i) = zero
706 funl(1)=ev(1,i)*yfac_unl*finter(ifunc(kk),ee(1,i),npf,tf,df)
707 funl(2)=ev(2,i)*yfac_unl*finter(ifunc(kk),ee(2,i),npf,tf,df)
708 funl(3)=ev(3,i)*yfac_unl*finter(ifunc(kk),ee(3,i),npf,tf,df)
709C
710 f(1)=ev(1,i)*yfac(1)*finter(ifunc(1),ee(1,i),npf,tf,df)
711 f(2)=ev(2,i)*yfac(1)*finter(ifunc(1),ee(2,i),npf,tf,df)
712 f(3)=ev(3,i)*yfac(1)*finter(ifunc(1),ee(3,i),npf,tf,df)
713C
714 IF(f(1)/=zero) dam(i) = max(dam(i),funl(1)/f(1))
715 IF(f(2)/=zero) dam(i) = max(dam(i),funl(2)/f(2))
716 IF(f(3)/=zero) dam(i) = max(dam(i),funl(3)/f(3))
717 IF(itype == -1 .OR. (itype == 2 .AND. rv(i) < one)) THEN
718 e_comp = min(zero,ee(1,i),ee(2,i),ee(3,i))
719 dam(i)= uvar(i,25)*dam(i)*abs(e_comp)
720 ELSE
721 e_tens = max(zero,ee(1,i),ee(2,i),ee(3,i))
722 dam(i)= uvar(i,22)*dam(i)*e_tens
723 ENDIF
724 !
725 fac = one
726 DO nn=1,20
727 fac = (-half)*fac
728 dd1 = ev(1,i)**fac - one
729 dd2 = ev(2,i)**fac - one
730 dd3 = ev(3,i)**fac - one
731 f(1) = f(1) +
732 . (dd1 + one)*yfac(1)*finter(ifunc(1),dd1,npf,tf,df)
733 f(2) = f(2) +
734 . (dd2 + one)*yfac(1)*finter(ifunc(1),dd2,npf,tf,df)
735 f(3) = f(3) +
736 . (dd3 + one)*yfac(1)*finter(ifunc(1),dd3,npf,tf,df)
737 ENDDO
738 invr = one/max(em20,rv(i))
739 p = invr*rbulk*(rv(i) - one)
740 t1(i) = invr*dam(i)*(two_third*f(1) - third*(f(2) + f(3))) + p
741 t2(i) = invr*dam(i)*(two_third*f(2) - third*(f(1) + f(3))) + p
742 t3(i) = invr*dam(i)*(two_third*f(3) - third*(f(1) + f(2))) + p
743 uvar(i,18) = dam(i)
744 ENDDO
745 ENDIF ! NLOAD
746 CASE(2) ! inunf_for = 2
747 ! Unloading with damage based on the energy
748 ! The unloading could be strain rate dependency
749 !
750 DO jj=1,ne_unload
751 i = indx_unl(jj)
752 dam(i) = zero
753 dam(i) = one - (uvar(i,16)/max(em20,uvar(i,17)))**shape
754 dam(i) = one - (one - hys)*dam(i)
755 ENDDO
756C
757 kk = nload + 1
758 DO jj=1,ne_unload
759 i = indx_unl(jj)
760 f(1)=ev(1,i)*yfac(1)*finter(ifunc(1),ee(1,i),npf,tf,df)
761 f(2)=ev(2,i)*yfac(1)*finter(ifunc(1),ee(2,i),npf,tf,df)
762 f(3)=ev(3,i)*yfac(1)*finter(ifunc(1),ee(3,i),npf,tf,df)
763 fac = one
764 DO nn=1,20
765 fac = (-half)*fac
766 dd1 = ev(1,i)**fac - one
767 dd2 = ev(2,i)**fac - one
768 dd3 = ev(3,i)**fac - one
769 f(1) = f(1) +
770 . (dd1 + one)*yfac(1)*finter(ifunc(1),dd1,npf,tf,df)
771 f(2) = f(2) +
772 . (dd2 + one)*yfac(1)*finter(ifunc(1),dd2,npf,tf,df)
773 f(3) = f(3) +
774 . (dd3 + one)*yfac(1)*finter(ifunc(1),dd3,npf,tf,df)
775 ENDDO
776 invr = one/max(em20,rv(i))
777 p = invr*rbulk*(rv(i) - one)
778 t1(i) = invr*dam(i)*(two_third*f(1) - third*(f(2) + f(3))) + p
779 t2(i) = invr*dam(i)*(two_third*f(2) - third*(f(1) + f(3))) + p
780 t3(i) = invr*dam(i)*(two_third*f(3) - third*(f(1) + f(2))) + p
781 uvar(i,18) = dam(i)
782 ENDDO
783C
784 CASE (3) !
785 ! Unloading with damage based on the energy
786 ! The unloading without dependency on strain rate
787 DO jj=1,ne_unload
788 i = indx_unl(jj)
789 dam(i) = zero
790 dam(i) = one - (uvar(i,16)/max(em20,uvar(i,17)))**shape
791 dam(i) = one - (one - hys)*dam(i)
792 ENDDO
793
794 DO jj=1,ne_unload
795 i = indx_unl(jj)
796C
797 IF(ee(1,i) >= zero) THEN
798 scale(1) = (uvar(i,22)* ee(1,i) + one)*yfac(1)
799 ELSE
800 scale(1) = (uvar(i,25)* ee(1,i) + one)*yfac(1)
801 ENDIF
802 IF(ee(2,i)>= zero) THEN
803 scale(2) = (uvar(i,22)* ee(2,i) + one)*yfac(1)
804 ELSE
805 scale(2) = (uvar(i,25)* ee(2,i) + one)*yfac(1)
806 ENDIF
807 IF(ee(3,i) >= zero) THEN
808 scale(3) = (uvar(i,22)* ee(3,i) + one)*yfac(1)
809 ELSE
810 scale(3) = (uvar(i,25)* ee(3,i) + one)*yfac(1)
811 ENDIF
812C
813 f(1) = ev(1,i)*scale(1)*finter(ifunc(1),ee(1,i),npf,tf,df1)
814 f(2) = ev(2,i)*scale(2)*finter(ifunc(1),ee(2,i),npf,tf,df2)
815 f(3) = ev(3,i)*scale(3)*finter(ifunc(1),ee(3,i),npf,tf,df3)
816
817 gmax(i) = max(gmax(i),scale(1)*df1,scale(2)*df2,
818 . scale(3)*df3 ) !
819 fac = one
820 DO nn=1,20
821 fac = fac*(-half)
822 dd1 = ev(1,i)**fac - one
823 dd2 = ev(2,i)**fac - one
824 dd3 = ev(3,i)**fac - one
825 f(1) = f(1) +
826 . (dd1 + one)*scale(1)*finter(ifunc(1),dd1,npf,tf,df)
827 f(2) = f(2) +
828 . (dd2 + one)*scale(2)*finter(ifunc(1),dd2,npf,tf,df)
829 f(3) = f(3) +
830 . (dd3 + one)*scale(3)*finter(ifunc(1),dd3,npf,tf,df)
831 ENDDO
832 invr = one/max(em20,rv(i))
833 p = invr*rbulk*(rv(i) - one)
834 t1(i) = invr*dam(i)*(two_third*f(1) - third*(f(2) + f(3))) + p
835 t2(i) = invr*dam(i)*(two_third*f(2) - third*(f(1) + f(3))) + p
836 t3(i) = invr*dam(i)*(two_third*f(3) - third*(f(1) + f(2))) + p
837 uvar(i,18) = dam(i)
838 ENDDO ! NE_UNLOAD
839 END SELECT
840 ENDIF ! NE_UNLOAD
841C
842C cauchy to glabale
843C
844 DO i=1,nel
845 signxx(i) = dirprv(1,1,i)*dirprv(1,1,i)*t1(i)
846 . + dirprv(1,2,i)*dirprv(1,2,i)*t2(i)
847 . + dirprv(1,3,i)*dirprv(1,3,i)*t3(i)
848
849 signyy(i) = dirprv(2,2,i)*dirprv(2,2,i)*t2(i)
850 . + dirprv(2,3,i)*dirprv(2,3,i)*t3(i)
851 . + dirprv(2,1,i)*dirprv(2,1,i)*t1(i)
852
853 signzz(i) = dirprv(3,3,i)*dirprv(3,3,i)*t3(i)
854 . + dirprv(3,1,i)*dirprv(3,1,i)*t1(i)
855 . + dirprv(3,2,i)*dirprv(3,2,i)*t2(i)
856
857 signxy(i) = dirprv(1,1,i)*dirprv(2,1,i)*t1(i)
858 . + dirprv(1,2,i)*dirprv(2,2,i)*t2(i)
859 . + dirprv(1,3,i)*dirprv(2,3,i)*t3(i)
860
861 signyz(i) = dirprv(2,2,i)*dirprv(3,2,i)*t2(i)
862 . + dirprv(2,3,i)*dirprv(3,3,i)*t3(i)
863 . + dirprv(2,1,i)*dirprv(3,1,i)*t1(i)
864
865 signzx(i) = dirprv(3,3,i)*dirprv(1,3,i)*t3(i)
866 . + dirprv(3,1,i)*dirprv(1,1,i)*t1(i)
867 . + dirprv(3,2,i)*dirprv(1,2,i)*t2(i)
868C
869* SET SOUND SPEED
870 soundsp(i)=sqrt((four_over_3*gs + rbulk)/rho(i))
871* SET VISCOSITY
872 viscmax(i) = zero
873 ENDDO
874 IF (impl_s > 0 .OR. ihet > 1) THEN
875 DO i=1,nel
876 et(i)= max(one,gmax(i)/gs)
877 ENDDO
878 ENDIF
879 RETURN
my_real
#define my_real
Definition cppsort.cpp:32
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
valpvecdp
subroutine valpvecdp(sig, val, vec, nel)
Definition matrix.F:533
valpvec
subroutine valpvec(sig, val, vec, nel)
Definition matrix.F:215