44
45
46
47 USE matparam_def_mod
49
50
51
52#include "implicit_f.inc"
53
54
55
56#include "mvsiz_p.inc"
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99#include "param_c.inc"
100
101 INTEGER ,INTENT(IN) :: NVARTMP,IMAT
102 INTEGER ,INTENT(INOUT) :: VARTMP(NEL,NVARTMP)
103 INTEGER NEL, NUPARAM, NUVAR,IPT,
104 . NGL(NEL),MAT(NEL),IPLA,IPM(NPROPMI,*)
106 . time,timestep,uparam(*),
107 . rho(nel),rho0(nel),volume(nel),eint(nel),
108 . depsxx(nel),depsyy(nel),depszz(nel),
109 . depsxy(nel),depsyz(nel),depszx(nel),
110 . epsxx(nel) ,epsyy(nel) ,epszz(nel) ,
111 . epsxy(nel) ,epsyz(nel) ,epszx(nel) ,
112 . sig0xx(nel),sig0yy(nel),sig0zz(nel),
113 . sig0xy(nel),sig0yz(nel),sig0zx(nel),
114 . pm(npropm,*)
115 TYPE(MATPARAM_STRUCT_) ,INTENT(INOUT) :: MATPARAM
116
117
118
120 . signxx(nel),signyy(nel),signzz
121 . signxy(nel),signyz(nel),signzx(nel),
122 . sigvxx(nel),sigvyy(nel),sigvzz(nel),
123 . sigvxy(nel),sigvyz(nel),sigvzx(nel),
124 . soundsp(nel),viscmax(nel)
125
126
127
129 . uvar(nel,nuvar), off(nel), pla(nel)
130
131
132
133 INTEGER NPF(*), MFUNC, KFUNC(MFUNC)
135 . tf(*),finter
136 EXTERNAL finter
137
138
139
140
141
142
143
144
145
146
147
148 INTEGER I,J,J1,J2,I1,I2,II,IADBUF,NC,NFUNC,ILOAD(MVSIZ),
149 . IFUNC(MVSIZ,100),NLOAD,,
150 . IE_CST(MVSIZ),ILOAD0,IDAMAGE,IFLAG,IPFUNC,NUMTABL
151 INTEGER :: NINDEX_EPSSM,NINDEX_EPST
152 INTEGER, DIMENSION(MVSIZ) :: INDEX_EPSSM,INDEX_EPST
154 . svm,inter_0,emax,epssmax,eps_max,dsig,expo,hys,de,df,mu
155 my_real yldelas(mvsiz),yldmin(mvsiz),yldmax(mvsiz),epss(mvsiz
156 . e0(mvsiz),g(mvsiz),nu(mvsiz)
157 . epst(mvsiz),aa(mvsiz),rate(mvsiz,2),
158 . yfac(mvsiz,2),sig0(mvsiz),eps0
159 INTEGER, DIMENSION(NEL) :: IPOS1
160 my_real,
DIMENSION(NEL) :: xvec1
161 INTEGER, DIMENSION(NEL,2) :: IPOS
162 my_real,
DIMENSION(NEL,2) :: xvec
163 my_real,
DIMENSION(MVSIZ) :: yld,dfc
164
165
166
167
168 nfunc = ipm(10,imat)
169 DO i= 1, nel
170 DO j=1,nfunc
171 ifunc(i,j)=ipm(10+j,imat)
172 ENDDO
173 ENDDO
174
175 iadbuf = ipm(7,imat)-1
176 nload = uparam(iadbuf+7)
177 nunload = uparam(iadbuf+8)
178 inter_0 = uparam(2*nfunc + 15)
179 epssmax = uparam(iadbuf+4)
180 DO i=1,nel
181 e0(i) = uparam(iadbuf+2)
182 aa(i) = uparam(iadbuf+3)
183 g(i) = uparam(iadbuf+5)
184 nu(i) = uparam(iadbuf+6)
185 idamage = uparam(iadbuf+ 2*nfunc + 9)
186 expo = uparam(iadbuf+ 2*nfunc + 10)
187 hys = uparam(iadbuf+ 2*nfunc + 11)
188 emax = uparam(iadbuf+ 2*nfunc + 12)
189 iflag = uparam(iadbuf+ 2*nfunc + 13)
190 yfacc = uparam(iadbuf+ 8 + 2*nfunc )
191 ENDDO
192 numtabl = matparam%NTABLE
193
194 DO i=1,nel
195
196
197
198 epst(i) = epsxx(i)**2+epsyy(i)**2 + epszz(i)**2 +
199 . half*(epsxy(i)**2+epsyz(i)**2+epszx(i)**2)
200 epst(i) = sqrt(epst(i))
201 epst(i) =
min(one, epst(i))
202
203
204
205 eps0(i) = uvar(i,1)
206 sig0(i) = uvar(i,2)
207 ENDDO
208
209
210
211 nindex_epst = 0
212 nindex_epssm = 0
213 DO i=1,nel
214
215 IF (epst(i) >= epssmax) THEN
216 yldelas(i)= inter_0
217 yldelas(i)= yldelas(i) + emax*(epst(i) - epssmax)
218 nindex_epssm = nindex_epssm + 1
219 index_epssm(nindex_epssm) = i
220 ELSE
221 nindex_epst = nindex_epst + 1
222 index_epst(nindex_epst) = i
223 ENDIF
224 ENDDO
225
226 IF (nindex_epst > 0) THEN
227 IF (matparam%TABLE(1)%NDIM == 1) THEN
228 DO ii=1,nindex_epst
229 i = index_epst(ii)
230 ipos1(i) = vartmp(i,1)
231 xvec1(i) = epst(i)
232 ENDDO
233
235
236 DO ii=1,nindex_epst
237 i = index_epst(ii)
238 yldelas(i) = yld(ii)
239 vartmp(i,1) = ipos1(ii)
240 ENDDO
241
242 ELSE
243 DO ii=1,nindex_epst
244 i = index_epst(ii)
245 ipos(ii,1) = vartmp(i,1)
246 ipos(ii,2) = 1
247 xvec(ii,1) = epst(i)
248 xvec(ii,2) = uparam(iadbuf + 9)
249 ENDDO
250
252
253 DO ii=1,nindex_epst
254 i = index_epst(ii)
255 yldelas(i) = yld(ii)
256 vartmp(i,1) = ipos(ii,1)
257 ENDDO
258 ENDIF
259
260 ENDIF
261
262
263
264
265 IF (matparam%TABLE(1)%NDIM == 1) THEN
266 DO i=1,nel
267 ipos1(i) = vartmp(i,1)
268 xvec1(i) = epst(i)
269 ENDDO
270
272
273 ELSE
274 DO i=1,nel
275 ipos(i,1) = vartmp(i,1)
276 ipos(i,2) = 1
277 xvec(i,1) = epst(i)
278 xvec(i,2) = uparam(iadbuf + 9)
279 ENDDO
280
282
283 END IF
284
285
286
287 IF (nunload > 0) THEN
288 IF (matparam%TABLE(2)%NDIM == 1) THEN
289 DO i=1,nel
290 ipos1(i) = 1
291 xvec1(i) = epst(i)
292 ENDDO
293
295
296 ELSE
297 DO i=1,nel
298 ipos(i,1) = 1
299 ipos(i,2) = 1
300 xvec(i,1) = epst(i)
301 xvec(i,2) = uparam(iadbuf + 9 + nload)
302 ENDDO
303
305
306 END IF
307 ELSE
308 yldmin(1:nel) = yldmax(1:nel)
309 END IF
310
311 IF(iflag > 0) THEN
312 DO i=1,nel
314 sig0xx(i) =
alpha*sig0xx(i)
315 sig0yy(i) =
alpha*sig0yy(i)
316 sig0zz(i) =
alpha*sig0zz(i)
317 sig0xy(i) =
alpha*sig0xy(i)
318 sig0zx(i) =
alpha*sig0zx(i)
319 sig0yz(i) =
alpha*sig0yz(i)
320 ENDDO
321 ENDIF
322
323
324 DO i = 1,nel
325 iload0 = uvar(i,5)
326 ie_cst(i)= 0
327 delta = epst(i) - uvar(i,4)
328 e = uvar(i,3)
329 uvar(i,9) = uvar(i,4)
330 IF(delta >= zero)THEN
331 yld(i) = yldmax(i)
332 iload(i) = 1
333 ELSEIF(delta < zero)THEN
334 yld(i) = yldmin(i)
335 iload(i) = -1
336 IF(idamage /= 0 )THEN
337 yld(i) = yldmax(i)
338 ENDIF
339 ENDIF
340
341 epss(i) = epst(i) - yld(i)/ e
342 epss(i) =
max(zero, epss(i))
343 de = aa(i)*(epss(i) - uvar(i,1))
344
345 IF(iload(i) == 1) THEN
346 e = e +
max(de, zero)
347 IF(iload0 == -1) e= uvar(i,3)
348 uvar(i,1) =
max(uvar(i,1), epss(i))
349 ELSE
350 e = e +
min(de ,zero)
351 IF(iload0 == 1) e= uvar(i,3)
352 uvar(i,1) =
min(epss(i),uvar(i,1))
353 ENDIF
356 uvar(i,3) = e
357
358 aa1(i) = e*(one-nu(i))/(one + nu(i))/(one - two*nu(i))
359 aa2(i) = aa1(i)*nu(i)/(one - nu(i))
360 g(i) =half*e/(one + nu(i))
361
362 signxx(i)= aa1(i)*depsxx(i) + aa2(i)*(depsyy(i) + depszz(i))
363 signyy(i)= aa1(i)*depsyy(i) + aa2(i)*(depsxx(i) + depszz(i))
364 signzz(i)= aa1(i)*depszz(i) + aa2(i)*(depsxx(i) + depsyy(i))
365 signxy(i)= g(i) *depsxy(i)
366 signyz(i)= g(i) *depsyz(i)
367 signzx(i)= g(i) *depszx(i)
368 dsig = (signxx(i)**2 + signyy(i)**2 + signzz(i)**2)
369 . + two*(signxy(i)**2 + signzx(i)**2 + signyz(i)**2)
370 dsig =sqrt(dsig)
371
372 signxx(i)=sig0xx(i) + aa1(i)*depsxx(i)
373 . + aa2(i)*(depsyy(i) + depszz(i))
374 signyy(i)=sig0yy(i) + aa1(i)*depsyy(i)
375 . + aa2(i)*(depsxx(i) + depszz(i))
376 signzz(i)=sig0zz(i) + aa1(i)*depszz(i)
377 . + aa2(i)*(depsxx(i) + depsyy(i))
378 signxy(i)=sig0xy
379 signyz(i)=sig0yz(i) + g(i) *depsyz(i)
380 signzx(i)=sig0zx(i) + g(i) *depszx(i)
381
382 svm2 = (signxx(i)**2 + signyy(i)**2 + signzz(i)**2)
383 . + two*(signxy(i)**2 + signzx(i)**2 + signyz(i)**2)
384 svm =sqrt(svm2)
385
386 soundsp(i) = sqrt(aa1(i)/rho0(i))
387 viscmax(i) = zero
388
389 IF(idamage == 0 ) THEN
390 IF(delta >= zero ) THEN
391 yld(i) =
min(yldmax(i),svm)
392 ELSEIF(delta < zero) THEN
393 yld(i) =
max(yldmin(i),svm)
394 yld(i) =
min(yld(i),uvar(i,6))
395 ie_cst(i) = 1
396 IF( dsig > yldmin(i) .AND. dsig > svm)THEN
397 yld(i) = yldmin(i)
398 ie_cst(i) = 0
399 ENDIF
400 ENDIF
401 ELSE
402 yld(i) = yldmax(i)
403 IF(delta > zero .AND. svm < yldmax(i
404 uvar(i,8) = uvar(i,8) +
405 . half*(yld(i) + uvar(i,6))*delta
406 uvar(i,8) =
max(zero, uvar(i,8))
407 uvar(i,2) =
max(uvar(i,2) , uvar(i,8))
408 ENDIF
409
410 ENDDO
411
412
413
414
415
416 IF(idamage == 0 ) THEN
417 DO i=1,nel
418 signxx(i)= aa1(i)*epsxx(i) + aa2(i)*(epsyy
419 signyy(i)= aa1(i)*epsyy(i) + aa2(i)*(epsxx(i) + epszz(i))
420 signzz(i)= aa1(i)*epszz(i) + aa2(i)*(epsxx(i) + epsyy(i))
421 signxy(i)= g(i) *epsxy(i)
422 signyz(i)= g(i) *epsyz(i)
423 signzx(i)= g(i) *epszx(i)
424
425 svm2 = (signxx(i)**2 + signyy(i)**2 + signzz(i)**2)
426 . + two*(signxy(i)**2 + signzx(i)**2 + signyz(i)**2)
427 svm =sqrt(svm2)
428 r = yld(i)/
max(em20,svm)
429 signxx(i)=signxx(i)*r
430 signyy(i)=signyy(i)*r
431 signzz(i)=signzz(i)*r
432 signxy(i)=signxy(i)*r
433 signyz(i)=signyz(i)*r
434 signzx(i)=signzx(i)*r
435
436 IF(ie_cst(i) == 1) THEN
437 iload0 = int(uvar(i,5))
438 IF(iload0 /= iload(i))THEN
439 iload(i) = iload0
440 uvar(i,1) = eps0(i)
441 ENDIF
442 ENDIF
443 uvar(i,4) = epst(i)
444 uvar(i,2) = svm*r
445 uvar(i,5) = iload(i)
446 uvar(i,6) = yld(i)
447 ENDDO
448 ELSEIF(idamage == 1) THEN
449 DO i=1,nel
450
451 r = one
452 signxx(i)= aa1(i)*epsxx(i) + aa2(i)*(epsyy(i) + epszz(i))
453 signyy(i)= aa1(i)*epsyy(i) + aa2(i)*(epsxx(i) + epszz(i))
454 signzz(i)= aa1(i)*epszz(i) + aa2(i)*(epsxx
455 signxy(i)= g(i) *epsxy(i)
456 signyz(i)= g(i) *epsyz(i)
457 signzx(i)= g(i) *epszx(i)
458
459 svm2 = (signxx(i)**2 + signyy(i)**2 + signzz(i)**2)
460 . + two*(signxy(i)**2 + signzx(i)**2 + signyz(i)**2)
461 svm =sqrt(svm2)
462 r = (yld(i)/
max(em20,svm))
463 signxx(i)=signxx(i)*r
464 signyy(i)=signyy(i)*r
465 signzz(i)=signzz(i)*r
466 signxy(i)=signxy(i)*r
467 signyz(i)=signyz(i)*r
468 signzx(i)=signzx(i)*r
469
470 IF(iload(i) == -1) THEN
471 r = (yldmin(i)/
max(em20,yldelas(i)))
472 p = third*(signxx(i) + signyy(i) + signzz(i))
473 signxx(i)=(signxx(i) - p)*r + p
474 signyy(i)=(signyy(i) - p)*r + p
475 signzz(i)=(signzz(i) - p)*r + p
476 signxy(i)=signxy(i)*r
477 signyz(i)=signyz(i)*r
478 signzx(i)=signzx(i)*r
479 ENDIF
480
481 svm2 = (signxx(i)**2 + signyy(i)**2 + signzz(i)**2)
482 . + two*(signxy(i)**2 + signzx(i)**2 + signyz(i)**2)
483 svm =sqrt(svm2)
484 uvar(i,4) = epst(i)
485 uvar(i,5) = iload(i)
486 uvar(i,6) = yld(i)
487 ENDDO
488 ELSEIF(idamage == 2 ) THEN
489 DO i=1,nel
490
491 signxx(i)= aa1(i)*epsxx(i) + aa2(i)*(epsyy(i) + epszz(i))
492 signyy(i)= aa1(i)*epsyy(i) + aa2(i)*(epsxx(i) + epszz(i))
493 signzz(i)= aa1(i)*epszz(i) + aa2(i)*(epsxx(i) + epsyy(i))
494 signxy(i)= g(i) *epsxy(i)
495 signyz(i)= g(i) *epsyz(i)
496 signzx(i)= g(i) *epszx(i)
497
498 svm2 = (signxx(i)**2 + signyy(i)**2 + signzz(i)**2)
499 . + two*(signxy(i)**2 + signzx(i)**2 + signyz(i)**2)
500 svm =sqrt(svm2)
501 r = yld(i)/
max(em20,svm)
502 signxx(i)=signxx(i)*r
503 signyy(i)=signyy(i)*r
504 signzz(i)=signzz(i)*r
505 signxy(i)=signxy(i)*r
506 signyz(i)=signyz(i)*r
507 signzx(i)=signzx(i)*r
508
509 IF(iload(i) == -1) THEN
510 r = yldmin(i)/
max(em20,yldelas(i))
511 signxx(i)=signxx(i)*r
512 signyy(i)=signyy(i)*r
513 signzz(i)=signzz(i)*r
514 signxy(i)=signxy(i)*r
515 signyz(i)=signyz(i)*r
516 signzx(i)=signzx(i)*r
517 ENDIF
518
519 svm2 = (signxx(i)**2 + signyy(i)**2 + signzz(i)**2)
520 . + two*(signxy(i)**2 + signzx(i)**2 + signyz(i)**2)
521 svm =sqrt(svm2)
522 uvar(i,4) = epst(i)
523 uvar(i,5) = iload(i)
524 uvar(i,6) = yld(i)
525 ENDDO
526
527 ELSEIF(idamage == 3) THEN
528 DO i=1,nel
529
530 signxx(i)= aa1(i)*epsxx(i) + aa2(i)*(epsyy(i) + epszz(i))
531 signyy(i)= aa1(i)*epsyy(i) + aa2(i)*(epsxx(i) + epszz(i))
532 signzz(i)= aa1(i)*epszz(i) + aa2(i)*(epsxx(i) + epsyy(i))
533 signxy(i)= g(i) *epsxy(i)
534 signyz(i)= g(i) *epsyz(i)
535 signzx(i)= g(i) *epszx(i)
536
537 svm2 = (signxx(i)**2 + signyy(i)**2 + signzz(i)**2)
538 . + two*(signxy(i)**2 + signzx(i)**2 + signyz(i)**2)
539 svm =sqrt(svm2)
540 r = (yld(i)/
max(em20,svm))
541 signxx(i)=signxx(i)*r
542 signyy(i)=signyy(i)*r
543 signzz(i)=signzz(i)*r
544 signxy(i)=signxy(i)*r
545 signyz(i)=signyz(i)*r
546 signzx(i)=signzx(i)*r
547
548 IF(iload(i) == -1) THEN
549 r = one - (uvar(i,8)/
max(em20,uvar(i,2)))**expo
550 r = one - (one - hys)*r
551 p = third*(signxx(i) + signyy(i) + signzz(i))
552 signxx(i)=(signxx(i) - p)*r + p
553 signyy(i)=(signyy(i) - p)*r + p
554 signzz(i)=(signzz(i) - p)*r + p
555 signxy(i)=signxy(i)*r
556 signyz(i)=signyz(i)*r
557 signzx(i)=signzx(i)*r
558 ENDIF
559
560 uvar(i,4) = epst(i)
561 uvar(i,5) = iload(i)
562 uvar(i,6) = yld(i)
563 ENDDO
564 ELSEIF(idamage == 4) THEN
565 DO i=1,nel
566
567 signxx(i)= aa1(i)*epsxx(i) + aa2(i)*(epsyy(i) + epszz(i))
568 signyy(i)= aa1(i)*epsyy(i) + aa2(i)*(epsxx(i) + epszz(i))
569 signzz(i)= aa1(i)*epszz(i) + aa2(i)*(epsxx(i) + epsyy(i))
570 signxy(i)= g(i) *epsxy(i)
571 signyz(i)= g(i) *epsyz(i)
572 signzx(i)= g(i) *epszx(i)
573
574 svm2 = (signxx(i)**2 + signyy(i)**2 + signzz(i)**2)
575 . + two*(signxy(i)**2 + signzx(i)**2 + signyz(i)**2)
576 svm =sqrt(svm2)
577 r = (yld(i)/
max(em20,svm))
578 signxx(i)=signxx(i)*r
579 signyy(i)=signyy(i)*r
580 signzz(i)=signzz(i)*r
581 signxy(i)=signxy(i)*r
582 signyz(i)=signyz(i)*r
583 signzx(i)=signzx(i)*r
584
585 IF(iload(i) == -1) THEN
586 r = one - (uvar(i,8)/
max(em20,uvar(i,2)))**expo
587 r = one - (one - hys)*r
588 signxx(i)=signxx(i)*r
589 signyy(i)=signyy(i)*r
590 signzz(i)=signzz(i)*r
591 signxy(i)=signxy(i)*r
592 signyz(i)=signyz(i)*r
593 signzx(i)=signzx(i)*r
594 ENDIF
595 uvar(i,4) = epst(i)
596 uvar(i,5) = iload(i)
597 uvar(i,6) = yld(i)
598 ENDDO
599 ENDIF
600
601 IF(iflag > 0 ) THEN
602 ipfunc = ifunc(1,nfunc)
603 DO i=1,nel
604 mu = 1 - rho(i)/rho0(i)
606 IF(mu > 0) THEN
607 alpha = yfacc*finter(ipfunc,mu,npf,tf,df)
609 ENDIF
610 signxx(i) =
alpha*signxx(i)
611 signyy(i) =
alpha*signyy(i)
612 signzz(i) =
alpha*signzz(i)
613 signxy(i) =
alpha*signxy(i)
614 signzx(i) =
alpha*signzx(i)
615 signyz(i) =
alpha*signyz(i)
617 ENDDO
618 ENDIF
619
620 RETURN
subroutine table_mat_vinterp(table, dimx, nel, ipos, xx, yy, dydx)