OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps80.F
Go to the documentation of this file.
1Copyright> OpenRadioss
2Copyright> Copyright (C) 1986-2025 Altair Engineering Inc.
3Copyright>
4Copyright> This program is free software: you can redistribute it and/or modify
5Copyright> it under the terms of the GNU Affero General Public License as published by
6Copyright> the Free Software Foundation, either version 3 of the License, or
7Copyright> (at your option) any later version.
8Copyright>
9Copyright> This program is distributed in the hope that it will be useful,
10Copyright> but WITHOUT ANY WARRANTY; without even the implied warranty of
11Copyright> MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12Copyright> GNU Affero General Public License for more details.
13Copyright>
14Copyright> You should have received a copy of the GNU Affero General Public License
15Copyright> along with this program. If not, see <https://www.gnu.org/licenses/>.
16Copyright>
17Copyright>
18Copyright> Commercial Alternative: Altair Radioss Software
19Copyright>
20Copyright> As an alternative to this open-source version, Altair also offers Altair Radioss
21Copyright> software under a commercial license. Contact Altair to discuss further if the
22Copyright> commercial version may interest you: https://www.altair.com/radioss/.
23!||====================================================================
24!|| sigeps80 ../engine/source/materials/mat/mat080/sigeps80.F
25!||--- called by ------------------------------------------------------
26!|| mulaw ../engine/source/materials/mat_share/mulaw.F90
27!||--- calls -----------------------------------------------------
28!|| finter ../engine/source/tools/curve/finter.F
29!|| kirkaldykinetics ../engine/source/materials/mat/mat080/kirkaldykinetics.F
30!|| phasekinetic2 ../engine/source/materials/mat/mat080/phasekinetic2.F
31!|| table_vinterp ../engine/source/tools/curve/table_tools.F
32!||--- uses -----------------------------------------------------
33!|| interface_table_mod ../engine/share/modules/table_mod.F
34!|| table_mod ../engine/share/modules/table_mod.f
35!||====================================================================
36 SUBROUTINE sigeps80(
37 1 NEL, NUPARAM, NUVAR, MFUNC,
38 2 KFUNC, NPF, TF, TIME,
39 3 TIMESTEP,UPARAM, RHO0, RHO,
40 4 VOL, EINT, EPSPXX, EPSPYY,
41 5 EPSPZZ, EPSPXY, EPSPYZ, EPSPZX,
42 6 DEPSXX, DEPSYY, DEPSZZ, DEPSXY,
43 7 DEPSYZ, DEPSZX, EPSXX, EPSYY,
44 8 EPSZZ, EPSXY, EPSYZ, EPSZX,
45 9 SIGOXX, SIGOYY, SIGOZZ, SIGOXY,
46 A SIGOYZ, SIGOZX, SIGNXX, SIGNYY,
47 B SIGNZZ, SIGNXY, SIGNYZ, SIGNZX,
48 C SOUNDSP, VISCMAX, UVAR, OFF,
49 D NGL, IPM, MAT, EPSP,
50 E YLD, PLA, TABLE, TEMP,
51 F NVARTMP, VARTMP, TRDEPSTH,EINTTH,
52 G JTHE ,IDT_THERM, THEACCFACT)
53C-----------------------------------------------
54 USE table_mod
55 USE interface_table_mod
56C-----------------------------------------------
57C I m p l i c i t T y p e s
58C-----------------------------------------------
59#include "implicit_f.inc"
60C---------+---------+---+---+--------------------------------------------
61C VAR | SIZE |TYP| RW| DEFINITION
62C---------+---------+---+---+--------------------------------------------
63C NEL | 1 | I | R | SIZE OF THE ELEMENT GROUP NEL F
64C NUPARAM | 1 | I | R | SIZE OF THE USER PARAMETER ARRAY
65C NUVAR | 1 | I | R | NUMBER OF USER ELEMENT VARIABLES
66C---------+---------+---+---+--------------------------------------------
67C NFUNC | 1 | I | R | NUMBER FUNCTION USED FOR THIS USER LAW
68C IFUNC | NFUNC | I | R | FUNCTION INDEX
69C NPF | * | I | R | FUNCTION ARRAY
70C NPT0 | 1 | I | R | NUMBER OF LAYERS OR INTEGRATION POINTS
71C IPT | 1 | I | R | LAYER OR INTEGRATION POINT NUMBER
72C IFLAG | * | I | R | GEOMETRICAL FLAGS
73C TF | * | F | R | FUNCTION ARRAY
74C---------+---------+---+---+--------------------------------------------
75C TIME | 1 | F | R | CURRENT TIME
76C TIMESTEP| 1 | F | R | CURRENT TIME STEP
77C UPARAM | NUPARAM | F | R | USER MATERIAL PARAMETER ARRAY
78C RHO0 | NEL | F | R | INITIAL DENSITY
79C AREA | NEL | F | R | AREA
80C EINT | 2*NEL | F | R | INTERNAL ENERGY(MEMBRANE,BENDING)
81C THKLY | NEL | F | R | LAYER THICKNESS
82C EPSPXX | NEL | F | R | STRAIN RATE XX
83C EPSPYY | NEL | F | R | STRAIN RATE YY
84C ... | | | |
85C DEPSXX | NEL | F | R | STRAIN INCREMENT XX
86C DEPSYY | NEL | F | R | STRAIN INCREMENT YY
87C ... | | | |
88C EPSXX | NEL | F | R | STRAIN XX
89C EPSYY | NEL | F | R | STRAIN YY
90C ... | | | |
91C SIGOXX | NEL | F | R | OLD ELASTO PLASTIC STRESS XX
92C SIGOYY | NEL | F | R | OLD ELASTO PLASTIC STRESS YY
93C ... | | | |
94C---------+---------+---+---+--------------------------------------------
95C SIGNXX | NEL | F | W | NEW ELASTO PLASTIC STRESS XX
96C SIGNYY | NEL | F | W | NEW ELASTO PLASTIC STRESS YY
97C ... | | | |
98C SIGVXX | NEL | F | W | VISCOUS STRESS XX
99C SIGVYY | NEL | F | W | VISCOUS STRESS YY
100C ... | | | |
101C SOUNDSP | NEL | F | W | SOUND SPEED (NEEDED FOR TIME STEP)
102C VISCMAX | NEL | F | W | MAXIMUM DAMPING MODULUS(NEEDED FOR TIME STEP)
103C---------+---------+---+---+--------------------------------------------
104C THK | NEL | F |R/W| THICKNESS
105C PLA | NEL | F |R/W| PLASTIC STRAIN
106C UVAR |NEL*NUVAR| F |R/W| USER ELEMENT VARIABLE ARRAY
107C OFF | NEL | F |R/W| DELETED ELEMENT FLAG (=1. ON, =0. OFF)
108C---------+---------+---+---+--------------------------------------------
109C C o m m o n B l o c k s
110C-----------------------------------------------
111#include "param_c.inc"
112#include "com01_c.inc"
113#include "scr18_c.inc"
114C-----------------------------------------------
115C I N P U T A r g u m e n t s
116C-----------------------------------------------
117C
118 INTEGER NEL, NUPARAM, NUVAR,NVARTMP, NPT0, IPT,
119 . NGL(NEL),MAT(NEL),IPM(NPROPMI,*),ITHK
120 my_real TIME,TIMESTEP,UPARAM(*),
121 . AREA(NEL),RHO0(NEL),EINT(NEL,2),
122 . THKLY(NEL),PLA(NEL),RHO(NEL),
123 . EPSPXX(NEL),EPSPYY(NEL),EPSPZZ(NEL),
124 . EPSPXY(NEL),EPSPYZ(NEL),EPSPZX(NEL),
125 . DEPSXX(NEL),DEPSYY(NEL),DEPSZZ(NEL),
126 . DEPSXY(NEL),DEPSYZ(NEL),DEPSZX(NEL),
127 . EPSXX(NEL) ,EPSYY(NEL) ,EPSZZ(NEL) ,
128 . EPSXY(NEL) ,EPSYZ(NEL) ,EPSZX(NEL) ,
129 . SIGOXX(NEL),SIGOYY(NEL),SIGOZZ(NEL),
130 . SIGOXY(NEL),SIGOYZ(NEL),SIGOZX(NEL),
131 . VOL(NEL) ,TEMP(NEL),
132 . DIE(NEL),COEF(NEL)
133 INTEGER, INTENT(IN) :: JTHE
134 INTEGER, INTENT(IN) :: IDT_THERM
135 my_real ,intent(in) :: THEACCFACT
136C-----------------------------------------------
137C O U T P U T A r g u m e n t s
138C-----------------------------------------------
139 my_real
140 . signxx(nel),signyy(nel),signzz(nel),
141 . signxy(nel),signyz(nel),signzx(nel),
142 . soundsp(nel),viscmax(nel)
143C-----------------------------------------------
144C I N P U T O U T P U T A r g u m e n t s
145C-----------------------------------------------
146 my_real
147 . uvar(nel,nuvar),epsp(nel),
148 . off(nel),thk(nel),yld(nel),eintth(nel)
149 INTEGER :: VARTMP(NEL,NVARTMP)
150
151 TYPE(ttable) TABLE(*)
152C-----------------------------------------------
153C VARIABLES FOR FUNCTION INTERPOLATION
154C-----------------------------------------------
155 INTEGER NPF(*), MFUNC, KFUNC(MFUNC)
156 my_real FINTER ,TF(*)
157 EXTERNAL FINTER
158c-----------------------------------------------
159C L o c a l V a r i a b l e s
160C-----------------------------------------------
161 INTEGER I,J,K,NRATE(NEL),II,J1,J2,N,NINDX,IADBUF,NINDXGLOB,NINDXLOC,
162 . nmax,index(nel),k2,iterk2,iter,niter,npfg(2),itable(5),
163 . efunc,ifunc(5),iterk,flag,ntab,heatflag,heatoption,
164 . niheat,nicool,flag_loc ,local,flag_tr_strain,flag_tr_kinetics
165
166 INTEGER, DIMENSION(NEL) :: INDXLOC,INDXGLOB
167 my_real
168 . cun ,cdeux,ctrois,deno,efac,
169 . alambda,blambda,ceps,peps,sol1,sol2,
170 . yscale(5),teta2,teta3,teta4,teta5,qr2,qr3,qr4,pp,
171 . alfa1,alfa2,t2,t3,t4,t5,t1,voli,xfac(5),rscale(5),
172 . ppxx,ppyy,ppzz,ppxy,ppyz,ppzx,esfim1,
173 . alpha,tref,ae1,ae3,bs,ms,gsize,nu,fcfer,fcper,fcbai,fg,
174 . e1po(nel),rhnew,eode,fgrain,kper,hfctn,kbain,xeq2,
175 . xeq4,lat1,lat2,hfp,hb,hm,vol0,tini,unitt,
176 . tgz(12),crit,test, heatflag1,dxrho,
177 . conv,huitcent,sspshell,sspsol,qptt,slope,dydxr,
178 . tau1, tau3,
179 . gfac_f,phi_f,psi_f,cr_f,cf,gfac_p,phi_p,psi_p,cr_p,cp,
180 . gfac_b,phi_b,psi_b,cr_b,cb,phi_m,psi_m,n_m,fgfer,fgper,fgbai,
181 . fac,dsvm,
182 . normxx,normyy,normzz,normxy,normyz,normzx,ddlam,dfdsig_dsig,dpla_dlam,ddep
183
184 my_real
185 . e(nel),rbulk(nel),shear(nel),g2(nel),tempel(nel),
186 . lamda(nel),dydxgz(nel),dydxe(nel),
187 . frac1(nel),frac2(nel),frac3(nel) ,frac4(nel),
188 . frac5(nel),totfrac(nel),totfracold(nel) ,soxx(nel) ,
189 . gz(nel) ,treo(nel) ,sigom(nel),soyy(nel),sozz(nel),
190 . dydx1(nel),dydx2(nel) ,dydx3(nel),
191 . dydx4(nel),dydx5(nel) ,yield(nel),
192 . y1(nel) ,y2(nel),y3(nel),y4(nel),y5(nel),
193 . strialxx(nel),strialyy(nel),strialzz(nel),
194 . strialxy(nel),strialyz(nel),strialzx(nel),svm(nel),
195 . eoxx(nel),eoyy(nel),eozz(nel),epsth(nel),
196 . eoxy(nel),eoyz(nel),eozx(nel),
197 . eodexx(nel),eodeyy(nel),eodezz(nel),
198 . eodexy(nel),eodezx(nel),eodeyz(nel),y1ini(nel),
199 . epsoxx(nel),epsoyy(nel),epsozz(nel),
200 . epsoxy(nel),epsoyz(nel),epsozx(nel),
201 . dpla(nel) ,snorm(nel),seff(nel),
202 . sigm(nel) ,sxx(nel),syy(nel),szz(nel),
203 . sxy(nel) ,syz(nel),szx(nel) ,
204 . deplxx(nel),deplyy(nel),deplzz(nel),deplxy(nel),deplyz(nel),deplzx(nel),
205 . dpla1(nel),dpla2(nel),dpla3(nel),dpla4(nel),dpla5(nel),
206 . pla1(nel) ,pla2(nel) ,pla3(nel) ,pla4(nel) ,pla5(nel),
207 . depselzz(nel),depsplzz(nel),depsim1(nel),depsi(nel),
208 . depsip1(nel),sigzim1(nel),sigi(nel),
209 . volold(nel),rh(nel),de12(nel),yldold(nel),
210 . depselozz(nel),depsplozz(nel),
211 . svmi(nel),trdepsth(nel),svmim1(nel),svmtr(nel),trdeps,
212 . depsthxx(nel),depsthyy(nel),depsthzz(nel),depsth(nel),
213 . x1(nel),x2(nel),x3(nel),x4(nel),x5(nel),hard(nel),
214 . dexx(nel),deyy(nel),dezz(nel),gold(nel),svmtest(nel),
215 . eplxx(nel),eplyy(nel),eplzz(nel) ,eplxy(nel),eplyz(nel),eplzx(nel),
216 . eeloxx(nel),eeloyy(nel),eeloyz(nel),eelozx(nel),
217 . eeloxy(nel),eelozz(nel),sigkk(nel),
218 . xgama(nel),tempmin(nel),vr(nel),h(nel),tempo(nel),
219 . zeq(nel), tau(nel),normdev(nel),deth12(nel),depstr(nel),
220 . rho_a(nel),rho_f(nel),rho_p(nel),rho_b(nel),rho_m(nel),rho_n(nel)
221 my_real, DIMENSION(NEL) ::
222 . a1, a2,fct,df,svmo,depsvol,hfct,a,b,c,dh,gsig,dgsig,esf,desf,p,svmt,
223 . lnf2,lnf3,lnf4,lnf5,faccs,stxx,styy,stzz,stxy,styz,stzx,logz,logzm1,
224 . dpxx,dpyy,dpzz,dpxy,dpyz,dpzx,phi2,phi3,phi4,phi5,psi2 ,psi3,psi4,psi5
225 my_real,
226 . DIMENSION(NEL,3) :: xx5,xx1,xx2,xx3,xx4
227
228 INTEGER,DIMENSION(NEL,3) :: IPOS1,IPOS2,IPOS3,IPOS4,IPOS5
229
230
231c-----------------------------------------------
232c USER VARIABLES INITIALIZATION
233c-----------------------------------------------
234C ! 0 < Ac1 < Ac3
235c UVAR 1= TEMPMIN
236c UVAR 2= FRAC1
237c UVAR 3= FRAC2
238c UVAR 4= FRAC3
239c UVAR 5= FRAC4
240c UVAR 6= FRAC5
241c UVAR 7= HARD
242c UVAR 8= TEMPEL
243c UVAR 9= YIELD
244c UVAR 10= XGAMA
245c UVAR 11= EPSXX
246c UVAR 12= EPSYY
247c UVAR 13= EPSZZ
248c UVAR 14= EPSXY
249c UVAR 15= SVM
250c UVAR 16= TIME
251c UVAR 17= PLA1
252c UVAR 18= PLA2
253c UVAR 19= PLA3
254c UVAR 20= PLA4
255c UVAR 21= PLA5
256c UVAR 22= TOTFRAC
257c UVAR 23= X2
258c UVAR 24= X3
259c UVAR 25= X4
260c UVAR 26= TIME
261c UVAR 27= TEMPEL
262c UVAR 29= EPLXX
263c UVAR 30= EPLYY
264c UVAR 31= EPLZZ
265c UVAR 32= EPLXY
266c UVAR 33= TEMPE
267c UVAR 34= TEMPE
268c UVAR 35= VOL
269c UVAR 36= SHEAR
270c UVAR 37= SIGNZZ
271c-----------------------------------------------
272
273 niter = 5
274 ntab = ipm(226,mat(1))
275
276 DO i=1,ntab
277 itable(i) = ipm(226+i,mat(1))
278 xfac(i) = uparam(58+i)
279 ENDDO
280 DO i=1,nel
281 e(i) = uparam(1)
282 ENDDO
283
284
285 DO j=1,5
286 yscale(j)=uparam(9+j)
287 ENDDO
288
289 nu = uparam(2)
290 efac = uparam(4)
291 ceps = uparam(15)
292 peps = uparam(16)
293 teta2 = uparam(17)
294 teta3 = uparam(18)
295 teta4 = uparam(19)
296 teta5 = uparam(20)
297 qr2 = uparam(21)
298 qr3 = uparam(22)
299 qr4 = uparam(23)
300 alpha = uparam(24)
301 tref = uparam(25)
302 ae1 = uparam(26)
303 ae3 = uparam(27)
304
305 bs = uparam(28)
306 ms = uparam(29)
307 gsize = uparam(30)
308
309 IF (idt_therm == 1) THEN ! ignore thermal expansion
310 alfa1 = zero
311 alfa2 = zero
312 ELSE
313 alfa1 = uparam(31)
314 alfa2 = uparam(32)
315 ENDIF
316 ! computed in starter
317 fcfer = uparam(33)
318 fcper = uparam(34)
319 fcbai = uparam(35)
320 fgrain = uparam(36)
321 kper = uparam(37)
322 kbain = uparam(38)
323 xeq2 = uparam(39)
324 lat1 = uparam(40)
325 lat2 = uparam(41)
326 hfp = uparam(42)
327 hb = uparam(43)
328 hm = uparam(44)
329 tini = uparam(45)
330 unitt = uparam(46)
331
332 gfac_f = uparam(75)
333 phi_f = uparam(76)
334 psi_f = uparam(77)
335 cr_f = uparam(78)
336
337 gfac_p = uparam(79)
338 phi_p = uparam(80)
339 psi_p = uparam(81)
340 cr_p = uparam(82)
341
342 gfac_b = uparam(83)
343 phi_b = uparam(84)
344 psi_b = uparam(85)
345 cr_b = uparam(86)
346
347 phi_m = uparam(84)
348 psi_m = uparam(85)
349 n_m = uparam(86)
350
351 fgfer = uparam(87)
352 fgper = uparam(88)
353 fgbai = uparam(89)
354C
355 cf = uparam(90)
356 cp = uparam(91)
357 cb = uparam(92)
358C
359 DO j=46+1,58
360 tgz(j-58+12)=uparam(j) !GZ FUNCTION
361 ENDDO
362
363 heatoption = uparam(58 + ntab + 1)
364 tau1 = uparam(58 + ntab + 2)
365 tau3 = uparam(58 + ntab + 3)
366 flag_loc = uparam(58 + ntab + 4)
367
368 flag_tr_strain = uparam(58 + ntab + 5)
369 flag_tr_kinetics = uparam(58 + ntab + 6)
370 rscale(1) = uparam(58 + ntab + 7)
371 rscale(2) = uparam(58 + ntab + 8)
372 rscale(3) = uparam(58 + ntab + 9)
373 rscale(4) = uparam(58 + ntab +10)
374 rscale(5) = uparam(58 + ntab +11)
375
376
377 heatflag = heatoption
378 IF(heatoption == 2) THEN
379 heatflag1 = finter(kfunc(2),time,npf,tf,slope)
380 heatflag = int(heatflag1)
381 ENDIF
382c
383 huitcent= eight*ep02
384 qptt=four+three*(em01+em02)
385 npfg(1)=1
386 npfg(2)=13
387 iterk =5
388 iterk2 =5
389 cun = third/(one-two*nu)
390 cdeux = half/(one+nu)
391 ctrois = nu/(one+nu)/(one-two*nu)
392
393 IF (isigi == 0) THEN
394 IF (time == zero)THEN
395 DO i=1,nel
396 uvar(i,35)=vol(i)
397 uvar(i,43)=rho(i)
398 IF(heatflag == 1)THEN
399 uvar(i,2) = zero ! fraction of austenite
400 uvar(i,3) = one ! fraction of ferrite
401 ELSE
402 uvar(i,2) = one ! fraction of austenite
403 ENDIF
404 uvar(i,8) = tini ! initial temperature
405 uvar(i,1) = tini
406
407 IF (jthe == -1) uvar(i,8) = temp(i)
408 xx1(i,1) = zero
409 xx1(i,2) = tini
410 xx1(i,3) = zero
411 ipos1(i,1) = 0
412 ipos1(i,2) = 0
413 ipos1(i,3) = 0
414 ENDDO
415!
416 CALL table_vinterp(table(itable(1)),nel,nel,ipos1,xx1,y1,dydx1)! initial yield
417!
418 uvar(1:nel,9) = y1(1:nel)
419 ENDIF
420 ENDIF
421!
422 IF (jthe /= 0 ) THEN
423 tempel(1:nel) = temp(1:nel)
424 ELSE
425 tempel(1:nel) = tini
426 ENDIF
427 DO i=1,nel
428 tempo(i) = uvar(i,8)
429 tempmin(i) = uvar(i,1) ! ASSURER PHASE CHANGE IF T DECREASE - avoid ossillations
430 ENDDO
431
432c-------------------------------------
433c compute temperature
434c-------------------------------------
435C-JTHE = -1 : thermal lagrangian, 0 : no thermal lagrange, 1 : thermal ale
436c IF(JTHE <= 0 ) THEN
437c DO I=1,NEL
438c IF (UVAR(I,8)<=EP03)THEN
439c CP(I)= 1.117*EP06/(1010.0-TEMPO(I))/(1010.0-TEMPO(I))+12622.
440c . /(1010.0 - TEMPO(I) )+ 0.3485*TEMPO(I)+ 355.6
441c ELSE
442c CP(I)= 1.225*EP08/((TEMPO(I)-990.)**4)+0.1381*TEMPO(I) + 585.7
443c ENDIF
444c TEMP = UPARAM(45)
445c VOL0 = VOL(I) * RHO0(I)
446c ! TEMPEL(I) = TINI + (EINT(I,1)+ EINT(I,2))/VOL(I)/(RHO0(I)*CP(I))
447c ENDDO
448c ENDIF
449c-------------------------------------
450 IF(kfunc(1) > zero)THEN
451 DO i=1,nel
452 e(i) = finter(kfunc(1),tempel(i),npf,tf,dydxe(i))
453 e(i)= e(i) * efac
454 ENDDO
455 ENDIF
456c-------------------------------------
457 DO i=1,nel
458 rbulk(i) = e(i)/three/(one - two*nu)
459 shear(i) = e(i)*half/(one+nu)
460 lamda(i) = e(i)*ctrois
461 g2(i) = two*shear(i)
462 a1(i) = e(i)*(one-nu) /((one + nu)*(one - two*nu))
463 a2(i) = a1(i)*nu/(one - nu)
464 ENDDO
465c-------------------------------------
466 DO i=1,nel
467 x2(i)= uvar(i,23)
468 x3(i)= uvar(i,24)
469 x4(i)= uvar(i,25)
470 x5(i)= uvar(i,44)
471 frac1(i)= uvar(i,2)
472 frac2(i)= uvar(i,3)
473 frac3(i)= uvar(i,4)
474 frac4(i)= uvar(i,5)
475 frac5(i)= uvar(i,6)
476 totfracold(i) = frac2(i)+frac3(i)+frac4(i)+frac5(i) ! Z sum of product phase
477 pla1(i)= uvar(i,17)
478 pla2(i)= uvar(i,18)
479 pla3(i)= uvar(i,19)
480 pla4(i)= uvar(i,20)
481 pla5(i)= uvar(i,21)
482 gold(i)= uvar(i,36)
483 xgama(i)= uvar(i,10)
484 dpla(i) = zero
485 dpla1(i)= zero
486 dpla2(i)= zero
487 dpla3(i)= zero
488 dpla4(i)= zero
489 dpla5(i)= zero
490 dpla(i) = zero
491 dpla1(i)= zero
492 dpla2(i)= zero
493 dpla3(i)= zero
494 dpla4(i)= zero
495 dpla5(i)= zero
496 phi2(i) = zero
497 psi2(i) = zero
498 phi3(i) = zero
499 psi3(i) = zero
500 phi4(i) = zero
501 psi4(i) = zero
502 phi5(i) = zero
503 psi5(i) = zero
504 ENDDO
505 IF (heatflag == 1) THEN
506 teta2 = zero
507 teta3 = zero
508 teta4 = zero
509 teta5 = zero
510 ENDIF
511c-------------------------------------
512c PHASE CHANGE calculation
513c-------------------------------------
514 nicool = 0
515 IF(flag_loc == 2) THEN ! global treatment - same behavior per part
516 IF (heatflag == 1) THEN ! heating activated AE1<AE3
517 DO i=1,nel
518 IF(tempel(i)>= uvar(i,8).AND.frac2(i)>zero)THEN
519 zeq(i) = (tempel(i)-ae1)/(ae3-ae1)
520 tau(i) = tau1 + zeq(i) * ( tau3 - tau1)
521 zeq(i) = min( one ,max( zero, zeq(i)))
522 tau(i) = max( tau3,min( tau1, tau(i)))
523 frac1(i) = uvar(i,2) + max(zero,(zeq(i) - uvar(i,2)) ) * timestep*theaccfact / tau(i)
524 frac2(i) = one - frac1(i) - frac3(i)- frac4(i)- frac5(i)
525 tempmin(i) = uvar(i,8)! heating activated AE1<AE3
526 ENDIF
527 ENDDO
528 ELSE ! HEATFLAG /= 1
529 DO i=1,nel
530 nicool = nicool + 1
531 index(nicool)=i
532 IF (tempel(i)<= 1073.0 .AND. uvar(i,26)==zero)THEN
533 uvar(i,26)=time
534 uvar(i,27)=tempel(i)
535 ENDIF
536 IF (tempel(i)<= 773.0 .AND. uvar(i,16)==zero)THEN
537 uvar(i,16)=time
538 uvar(i,33)=tempel(i)
539 ENDIF
540 ENDDO
541
542 IF (flag_tr_kinetics == 2) THEN
543 CALL phasekinetic2(nel,time,tempel,tempo,tempmin,ae1,ae3,bs,ms,fcfer,fcper,
544 . fcbai,fgrain,frac1,frac2,frac3,frac4,frac5,x2,x3,x4,x5,
545 . gfac_f,phi_f,psi_f,cr_f,cf,gfac_p,phi_p,psi_p,cr_p,cp,
546 . gfac_b,phi_b,psi_b,cr_b,cb,phi_m,psi_m,n_m,fgfer,fgper,fgbai,
547 . qr2,qr3,qr4,kper,kbain,alpha,xeq2,xeq4,xgama,totfracold,timestep,nicool,
548 . index,theaccfact)
549 ELSE
550 CALL kirkaldykinetics(nel,time,tempel,tempmin,ae1,ae3,bs,ms,fcfer,fcper,
551 . fcbai,fgrain,frac1,frac2,frac3,frac4,frac5,x2,x3,x4,x5,qr2,
552 . qr3,qr4,kper,kbain,alpha,xeq2,xeq4,xgama,totfracold,timestep,nicool,
553 . index,theaccfact)
554 ENDIF
555 ENDIF
556 ELSE !FLAG_LOC /= 2 local
557
558 DO i=1,nel
559 IF(tempel(i)>= uvar(i,8).AND.frac2(i)>zero)THEN
560 zeq(i) = (tempel(i)-ae1)/(ae3-ae1)
561 tau(i) = tau1 + zeq(i) * ( tau3 - tau1)
562 zeq(i) = min( one ,max( zero, zeq(i)))
563 tau(i) = max( tau3,min( tau1, tau(i)))
564 frac1(i) = uvar(i,2) + max(zero,(zeq(i) - uvar(i,2)) ) * timestep*theaccfact / tau(i)
565 frac2(i) = one - frac1(i) - frac3(i)- frac4(i)- frac5(i)
566 tempmin(i) = uvar(i,8)! heating activated AE1<AE3
567 ELSE !cooling
568 nicool = nicool+1
569 index(nicool)=i
570 IF (tempel(i)<= 1073.0 .AND. uvar(i,26)==zero)THEN
571 uvar(i,26)=time
572 uvar(i,27)=tempel(i)
573 ENDIF
574
575 IF (tempel(i)<= 773.0 .AND. uvar(i,16)==zero)THEN
576 uvar(i,16)=time
577 uvar(i,33)=tempel(i)
578 ENDIF
579 ENDIF
580 ENDDO
581 IF (nicool > 0) THEN
582 IF (flag_tr_kinetics == 2) THEN
583 CALL phasekinetic2(nel,time,tempel,tempo,tempmin,ae1,ae3,bs,ms,fcfer,fcper,
584 . fcbai,fgrain,frac1,frac2,frac3,frac4,frac5,x2,x3,x4,x5,
585 . gfac_f,phi_f,psi_f,cr_f,cf,gfac_p,phi_p,psi_p,cr_p,cp,
586 . gfac_b,phi_b,psi_b,cr_b,cb,phi_m,psi_m,n_m,fgfer,fgper,fgbai,
587 . qr2,qr3,qr4,kper,kbain,alpha,xeq2,xeq4,xgama,totfracold,timestep,nicool,
588 . index ,theaccfact)
589 ELSE
590 CALL kirkaldykinetics(nel,time,tempel,tempmin,ae1,ae3,bs,ms,fcfer,fcper,
591 . fcbai,fgrain,frac1,frac2,frac3,frac4,frac5,x2,x3,x4,x5,qr2,
592 . qr3,qr4,kper,kbain,alpha,xeq2,xeq4,xgama,totfracold,timestep,nicool,
593 . index ,theaccfact)
594 ENDIF
595 ENDIF
596 endif!FLAG_LOC /= 2
597c-----------------------------------------------
598 DO i=1,nel
599 IF(tempo(i)<=tempmin(i))tempmin(i)=tempo(i)
600 uvar(i,1) = tempmin(i)
601 ENDDO
602c-----------------------------------------------
603c-------------------------------------
604c interpolation of yield for each phase
605c-------------------------------------
606 DO i=1,nel
607
608 ipos1(i,1) = vartmp(i,1)
609 ipos1(i,2) = vartmp(i,2)
610 ipos1(i,3) = vartmp(i,3)
611C
612 ipos2(i,1) = vartmp(i,4)
613 ipos2(i,2) = vartmp(i,5)
614 ipos2(i,3) = vartmp(i,6)
615
616 ipos3(i,1) = vartmp(i,7)
617 ipos3(i,2) = vartmp(i,8)
618 ipos3(i,3) = vartmp(i,9)
619C
620 ipos4(i,1) = vartmp(i,10)
621 ipos4(i,2) = vartmp(i,11)
622 ipos4(i,3) = vartmp(i,12)
623C
624 ipos5(i,1) = vartmp(i,13)
625 ipos5(i,2) = vartmp(i,14)
626 ipos5(i,3) = vartmp(i,15)
627
628C
629 xx1(i,1)=pla1(i)
630 xx1(i,2)=tempel(i)
631 xx1(i,3)=epsp(i)*xfac(1)
632C
633 xx2(i,1)=pla2(i)
634 xx2(i,2)=tempel(i)
635 xx2(i,3)=epsp(i)*xfac(2)
636C
637 xx3(i,1)=pla3(i)
638 xx3(i,2)=tempel(i)
639 xx3(i,3)=epsp(i)*xfac(3)
640C
641 xx4(i,1)=pla4(i)
642 xx4(i,2)=tempel(i)
643 xx4(i,3)=epsp(i)*xfac(4)
644C
645 xx5(i,1)=pla5(i)
646 xx5(i,2)=tempel(i)
647 xx5(i,3)=epsp(i)*xfac(5)
648 END DO
649C
650 CALL table_vinterp(table(itable(1)),nel,nel,ipos1,xx1,y1,dydx1)
651 CALL table_vinterp(table(itable(2)),nel,nel,ipos2,xx2,y2,dydx2)
652 CALL table_vinterp(table(itable(3)),nel,nel,ipos3,xx3,y3,dydx3)
653 CALL table_vinterp(table(itable(4)),nel,nel,ipos4,xx4,y4,dydx4)
654 CALL table_vinterp(table(itable(5)),nel,nel,ipos5,xx5,y5,dydx5)
655
656 DO i=1,nel
657
658 y1(i)=y1(i)*yscale(1)
659 y2(i)=y2(i)*yscale(2)
660 y3(i)=y3(i)*yscale(3)
661 y4(i)=y4(i)*yscale(4)
662 y5(i)=y5(i)*yscale(5)
663
664 dydx1(i)=dydx1(i)*yscale(1)
665 dydx2(i)=dydx2(i)*yscale(2)
666 dydx3(i)=dydx3(i)*yscale(3)
667 dydx4(i)=dydx4(i)*yscale(4)
668 dydx5(i)=dydx5(i)*yscale(5)
669 yield(i)=uvar(i,9)
670 yld(i)=max(em20,y1(i)*frac1(i)+y2(i)*frac2(i)+
671 . y3(i)*frac3(i)+y4(i)*frac4(i)+y5(i)*frac5(i))
672
673
674 y1ini(i) = max(em20,y1(i))
675 vartmp(i,1) = ipos1(i,1)
676 vartmp(i,2) = ipos1(i,2)
677 vartmp(i,3) = ipos1(i,3)
678C
679 vartmp(i,4) = ipos2(i,1)
680 vartmp(i,5) = ipos2(i,2)
681 vartmp(i,6) = ipos2(i,3)
682
683 vartmp(i,7) = ipos3(i,1)
684 vartmp(i,8) = ipos3(i,2)
685 vartmp(i,9) = ipos3(i,3)
686C
687 vartmp(i,10) = ipos4(i,1)
688 vartmp(i,11) = ipos4(i,2)
689 vartmp(i,12) = ipos4(i,3)
690C
691 vartmp(i,13) = ipos5(i,1)
692 vartmp(i,14) = ipos5(i,2)
693 vartmp(i,15) = ipos5(i,3)
694 ENDDO
695 faccs(1:nel) = one
696 ! analytic strain rate dependency
697 IF(ceps/=zero.AND.peps/=zero)THEN
698 DO i=1,nel
699 IF(epsp(i)>ceps)THEN
700 faccs(i)=one+ (epsp(i)/ceps)**(one/peps)
701 yld(i)= yld(i)*faccs(i)
702 ENDIF
703 ENDDO
704 ENDIF
705c---------------------------------------------- ------------------
706c----------------------------------------------------------------
707c compute thermal strain and transformation strain
708c---------------------------------------------- ------------------
709c----------------------------------------------------------------
710 DO i=1,nel
711 totfrac(i) = frac2(i)+frac3(i)+frac4(i)+frac5(i)
712 depsth(i) = (frac1(i)*alfa1+totfrac(i)*alfa2)*(tempel(i)-uvar(i,8))
713 trdepsth(i) = three*depsth(i)
714 ENDDO
715
716 IF(flag_tr_strain == 2) THEN
717 ! dependent directly on density
718 ! compute density for eachphase depending on temperature
719 ! interpole aust -f -p- b- m
720 DO i=1,nel
721 rho_a(i) = finter(kfunc(3),tempel(i),npf,tf,dydxr)*rscale(1)
722 rho_f(i) = finter(kfunc(4),tempel(i),npf,tf,dydxr)*rscale(2)
723 rho_p(i) = finter(kfunc(5),tempel(i),npf,tf,dydxr)*rscale(3)
724 rho_b(i) = finter(kfunc(6),tempel(i),npf,tf,dydxr)*rscale(4)
725 rho_m(i) = finter(kfunc(7),tempel(i),npf,tf,dydxr)*rscale(5)
726 dxrho = (frac1(i) - uvar(i,2)) * rho_a(i) +
727 . (frac2(i) - uvar(i,3)) * rho_f(i) +
728 . (frac3(i) - uvar(i,4)) * rho_p(i) +
729 . (frac4(i) - uvar(i,5)) * rho_b(i) +
730 . (frac5(i) - uvar(i,6)) * rho_m(i)
731 rho_n(i) = frac1(i)*rho_a(i) + frac2(i)*rho_f(i) + frac3(i)*rho_p(i) + frac4(i)*rho_b(i) + frac5(i)*rho_m(i)
732 depstr(i) = zero
733 IF(totfrac(i) > zero.AND.totfrac(i)< one)
734 . depstr(i) = third*dxrho/(rho0(i) + rho_n(i)-uvar(i,43))
735 uvar(i,43)= rho_n(i)
736 ENDDO
737 ELSE
738 DO i=1,nel
739 deth12(i)=qptt*em03
740 IF (tempel(i)<ms ) deth12(i)=six*em03
741 depstr(i) = zero
742 IF(totfrac(i) > zero.AND.totfrac(i)< one)
743 . depstr(i) = deth12(i) * (totfrac(i)-totfracold(i))!*LOG(TOTFRAC(I))/ (ONE - TOTFRAC(I))
744 ENDDO
745 ENDIF
746 DO i=1,nel
747 depsxx(i) = depsxx(i) - depsth(i) - depstr(i)
748 depsyy(i) = depsyy(i) - depsth(i) - depstr(i)
749 depszz(i) = depszz(i) - depsth(i) - depstr(i)
750 !UVAR(I,42) = UVAR(I,42) + DEPSTH(I)+ DEPSTR(I)
751
752 ENDDO
753 IF (timestep /=zero) THEN
754 DO i=1,nel
755 epspxx(i) = depsxx(i) /timestep
756 epspyy(i) = depsyy(i) /timestep
757 epspzz(i) = depszz(i) /timestep
758 ENDDO !DO I=1,NEL
759 ENDIF
760 DO i=1,nel
761 gz(i) = finter(1, totfrac(i),npfg,tgz,dydxgz(i))! used to compute plastic strain rate
762 ENDDO
763c----------------------------------------------
764 DO i=1,nel
765 IF (off(i)== one )THEN
766 lnf2(i) = zero
767 lnf3(i) = zero
768 lnf4(i) = zero
769 lnf5(i) = zero !used in plast strain increment local yield
770 IF(frac2(i)>uvar(i,3).AND. uvar(i,3) > zero) lnf2(i) = log(frac2(i)/uvar(i,3) )
771 IF(frac3(i)>uvar(i,4).AND. uvar(i,4) > zero) lnf3(i) = log(frac3(i)/uvar(i,4) )
772 IF(frac4(i)>uvar(i,5).AND. uvar(i,5) > zero) lnf4(i) = log(frac4(i)/uvar(i,5) )
773 IF(frac5(i)>uvar(i,6).AND. uvar(i,6) > zero) lnf5(i) = log(frac5(i)/uvar(i,6) )
774 ENDIF
775 ENDDO !DO I=1,NEL
776c----------------------------------------------
777c================================================
778c Mechanical calculation -+ DEPSTH(I)
779c================================================
780 DO i=1,nel
781 ! Computation of the trial stress tensor
782 signxx(i) = sigoxx(i) + a1(i)*depsxx(i)+ a2(i)*(depsyy(i) + depszz(i))
783 signyy(i) = sigoyy(i) + a1(i)*depsyy(i)+ a2(i)*(depsxx(i) + depszz(i))
784 signzz(i) = sigozz(i) + a1(i)*depszz(i)+ a2(i)*(depsxx(i) + depsyy(i))
785 signxy(i) = sigoxy(i) + shear(i)*depsxy(i)
786 signyz(i) = sigoyz(i) + shear(i)*depsyz(i)
787 signzx(i) = sigozx(i) + shear(i)*depszx(i)
788 ! Computation of the pressure of the trial stress tensor
789 p(i) = -third*(signxx(i) + signyy(i) + signzz(i))
790 ! Computation of the Von Mises equivalent stress of the trial stress tensor
791 sxx(i) = signxx(i) + p(i)
792 syy(i) = signyy(i) + p(i)
793 szz(i) = signzz(i) + p(i)
794 sxy(i) = signxy(i)
795 syz(i) = signyz(i)
796 szx(i) = signzx(i)
797 svm(i) = half*(sxx(i)**2 + syy(i)**2 + szz(i)**2)
798 . + (sxy(i)**2 + szx(i)**2 + syz(i)**2)
799 svm(i) = sqrt(three*svm(i))
800
801 normdev(i) = sqrt(sxx(i)* sxx(i)
802 . + syy(i)* syy(i)
803 . + szz(i)* szz(i) ) /e(i)
804
805
806 sigom(i) = -(sigoxx(i)+sigoyy(i)+sigozz(i)) * third
807 soxx(i) = sigoxx(i)+sigom(i)
808 soyy(i) = sigoyy(i)+sigom(i)
809 sozz(i) = sigozz(i)+sigom(i)
810 svmo(i) = sqrt(three_half*(soxx(i)*soxx(i)
811 . + soyy(i)* soyy(i)
812 . + sozz(i)* sozz(i)
813 . +two*(sigoxy(i)*sigoxy(i)
814 . +sigoyz(i)*sigoyz(i)
815 . +sigozx(i)*sigozx(i)) ) )
816 ENDDO
817c----------------------------------------------
818c----------------------------------------------
819 !-------------------------------------------------------------------------
820 ! check if local or global yield -in both cases plastic deformation occurs
821 !-------------------------------------------------------------------------
822 nindxglob = 0
823 nindxloc = 0
824 DO i=1,nel
825 IF ( svm(i) < yld(i) .AND. off(i) == one
826 . .AND.totfrac(i)>totfracold(i).AND. normdev(i)> em10
827 . .AND.totfrac(i)< one .AND. svm(i)>svmo(i)) THEN !.AND.HEATFLAG==0
828 ! LOCAL YIELD ALGORITHM - ONLY OCCURS IF DEVIATORIC STRESSES ARE APPLIED
829 !----------------------
830 nindxloc = nindxloc + 1
831 indxloc(nindxloc) = i
832 ENDIF
833 ENDDO
834 DO i=1,nel
835 IF (svm(i) > yld(i) .AND. off(i) == one) THEN
836 ! GLOBAL YIELD ALGORITHM
837 !----------------------
838 nindxglob = nindxglob + 1
839 indxglob(nindxglob) = i
840 ENDIF
841 ENDDO
842
843 !========================================================================================
844 ! GLOBAL YIELD ALGORITHM
845 !========================================================================================
846 IF (nindxglob > 0) THEN
847 DO ii = 1,nindxglob
848 i = indxglob(ii)
849 IF (uvar(i,3)/=zero)THEN
850 psi2(i) = max(zero,(lnf2(i)*(teta2*pla1(i)-pla2(i)))/(one+lnf2(i)) )
851 phi2(i) =(one+teta2*lnf2(i))/(one+lnf2(i))
852 ENDIF
853 IF (uvar(i,4)/=zero)THEN
854 psi3(i) = max(zero,(lnf3(i)*(teta3*pla1(i)-pla3(i)))/(one+lnf3(i)) )
855 phi3(i) =(one+teta3*lnf3(i))/(one+lnf3(i))
856 ENDIF
857 IF (uvar(i,5)/=zero)THEN
858 psi4(i) = max(zero, (lnf4(i)*(teta4*pla1(i)-pla4(i)))/(one+lnf4(i)) )
859 phi4(i) =(one+teta4*lnf4(i))/(one+lnf4(i))
860 ENDIF
861 IF (uvar(i,6)/=zero)THEN
862 psi5(i) = max(zero,(lnf5(i)*(teta5*pla1(i)-pla5(i)))/(one+lnf5(i)) )
863 phi5(i) = (one+teta5*lnf5(i))/(one+lnf5(i))
864 ENDIF
865 ENDDO !II = 1,NINDXGLOB
866 DO iter = 1,5
867 DO ii = 1,nindxglob
868 i = indxglob(ii)
869 fct(i) = svm(i) - yld(i)
870 !df/dsig
871 normxx = three_half * sxx(i) /svm(i) ! DF/DSIG
872 normyy = three_half * syy(i) /svm(i)
873 normzz = three_half * szz(i) /svm(i)
874 normxy = two * three_half * signxy(i) /svm(i)
875 normyz = two * three_half * signyz(i) /svm(i)
876 normzx = two * three_half * signzx(i) /svm(i)
877 ! df/ddlam = DF/DSIG * DSIG/DDLAM
878 df(i) = normxx * (a1(i)*normxx + a2(i)*(normyy+normzz))
879 . + normyy * (a1(i)*normyy + a2(i)*(normxx+normzz))
880 . + normzz * (a1(i)*normzz + a2(i)*(normxx+normyy))
881 . + normxy * normxy * shear(i)
882 . + normyz * normyz * shear(i)
883 . + normzx * normzx * shear(i)
884
885 df(i) = sign(max(abs(df(i)),em20) ,df(i))
886 ddlam = fct(i)/df(i)
887
888 !compute derivative of effective plastic strain vs dlam
889 dfdsig_dsig = signxx(i) * normxx + signyy(i) * normyy + signzz(i) * normzz
890 . + signxy(i) * normxy + signyz(i) * normyz + signzx(i) * normzx
891 dpla_dlam = dfdsig_dsig / yld(i)
892
893
894 ! Plastic strains tensor update
895 dpxx(i) = ddlam * normxx
896 dpyy(i) = ddlam * normyy
897 dpzz(i) = ddlam * normzz
898 dpxy(i) = ddlam * normxy
899 dpyz(i) = ddlam * normyz
900 dpzx(i) = ddlam * normzx
901 ! Elasto-plastic stresses update
902 signxx(i) = signxx(i) - (a1(i)*dpxx(i) + a2(i)*(dpyy(i) + dpzz(i)))
903 signyy(i) = signyy(i) - (a1(i)*dpyy(i) + a2(i)*(dpxx(i) + dpzz(i)))
904 signzz(i) = signzz(i) - (a1(i)*dpzz(i) + a2(i)*(dpxx(i) + dpyy(i)))
905 signxy(i) = signxy(i) - shear(i)*dpxy(i)
906 signyz(i) = signyz(i) - shear(i)*dpyz(i)
907 signzx(i) = signzx(i) - shear(i)*dpzx(i)
908 p(i) = -third*(signxx(i) + signyy(i) + signzz(i))
909 ! Computation of the Von Mises equivalent stress of the stress tensor
910 sxx(i) = signxx(i) + p(i)
911 syy(i) = signyy(i) + p(i)
912 szz(i) = signzz(i) + p(i)
913 sxy(i) = signxy(i)
914 syz(i) = signyz(i)
915 szx(i) = signzx(i)
916 svm(i) = half*(sxx(i)**2 + syy(i)**2 + szz(i)**2)
917 . + (sxy(i)**2 + szx(i)**2 + syz(i)**2)
918 svm(i) = sqrt(three*svm(i))
919
920 ddep = ddlam * dpla_dlam !DDEP IN THE ITERATIONS (DPLA_K+1 = DPLA_K + DDEP)
921 dpla(i) = max(zero, dpla(i) + ddep )
922
923 ! DPLA AND PLA OF EACH PHASE
924 dpla1(i)= dpla1(i) + ddep !DPLA(I)
925 IF (uvar(i,3)>zero)THEN
926 dpla2(i)= dpla1(i) *phi2(i) + psi2(i)
927 ENDIF
928 IF (uvar(i,4)>zero)THEN
929 dpla3(i)= dpla1(i) *phi3(i) + psi3(i)
930 ENDIF
931 IF (uvar(i,5)>zero)THEN
932 dpla4(i)= dpla1(i) *phi4(i) + psi4(i)
933 ENDIF
934 IF (uvar(i,6)>zero)THEN
935 dpla5(i)= dpla1(i) *phi5(i) + psi5(i)
936 ENDIF
937 !-------------------------------------
938 !UPDATE yield for each phase
939 !-------------------------------------
940 ipos1(i,1) = vartmp(i,1)
941 ipos1(i,2) = vartmp(i,2)
942 ipos1(i,3) = vartmp(i,3)
943C
944 ipos2(i,1) = vartmp(i,4)
945 ipos2(i,2) = vartmp(i,5)
946 ipos2(i,3) = vartmp(i,6)
947
948 ipos3(i,1) = vartmp(i,7)
949 ipos3(i,2) = vartmp(i,8)
950 ipos3(i,3) = vartmp(i,9)
951C
952 ipos4(i,1) = vartmp(i,10)
953 ipos4(i,2) = vartmp(i,11)
954 ipos4(i,3) = vartmp(i,12)
955C
956 ipos5(i,1) = vartmp(i,13)
957 ipos5(i,2) = vartmp(i,14)
958 ipos5(i,3) = vartmp(i,15)
959C
960 xx1(i,1)=pla1(i) + dpla1(i)
961 xx2(i,1)=pla2(i) + dpla2(i)
962 xx3(i,1)=pla3(i) + dpla3(i)
963 xx4(i,1)=pla4(i) + dpla4(i)
964 xx5(i,1)=pla5(i) + dpla5(i)
965 END DO
966C
967 CALL table_vinterp(table(itable(1)),nel,nel,ipos1,xx1,y1,dydx1)
968 CALL table_vinterp(table(itable(2)),nel,nel,ipos2,xx2,y2,dydx2)
969 CALL table_vinterp(table(itable(3)),nel,nel,ipos3,xx3,y3,dydx3)
970 CALL table_vinterp(table(itable(4)),nel,nel,ipos4,xx4,y4,dydx4)
971 CALL table_vinterp(table(itable(5)),nel,nel,ipos5,xx5,y5,dydx5)
972 DO ii = 1,nindxglob
973 i = indxglob(ii)
974 y1(i)=y1(i)*yscale(1)
975 y2(i)=y2(i)*yscale(2)
976 y3(i)=y3(i)*yscale(3)
977 y4(i)=y4(i)*yscale(4)
978 y5(i)=y5(i)*yscale(5)
979 dydx1(i)=dydx1(i)*yscale(1)
980 dydx2(i)=dydx2(i)*yscale(2)
981 dydx3(i)=dydx3(i)*yscale(3)
982 dydx4(i)=dydx4(i)*yscale(4)
983 dydx5(i)=dydx5(i)*yscale(5)
984
985 !new yield updated in newton iterations
986 yld(i)=max(em20, y1(i)*frac1(i)+
987 . y2(i)*frac2(i)+
988 . y3(i)*frac3(i)+
989 . y4(i)*frac4(i)+
990 . y5(i)*frac5(i))
991 yld(i)= yld(i)*faccs(i)
992
993 vartmp(i,1) = ipos1(i,1)
994 vartmp(i,2) = ipos1(i,2)
995 vartmp(i,3) = ipos1(i,3)
996C
997 vartmp(i,4) = ipos2(i,1)
998 vartmp(i,5) = ipos2(i,2)
999 vartmp(i,6) = ipos2(i,3)
1000
1001 vartmp(i,7) = ipos3(i,1)
1002 vartmp(i,8) = ipos3(i,2)
1003 vartmp(i,9) = ipos3(i,3)
1004C
1005 vartmp(i,10) = ipos4(i,1)
1006 vartmp(i,11) = ipos4(i,2)
1007 vartmp(i,12) = ipos4(i,3)
1008C
1009 vartmp(i,13) = ipos5(i,1)
1010 vartmp(i,14) = ipos5(i,2)
1011 vartmp(i,15) = ipos5(i,3)
1012 ENDDO ! indexed elements
1013 enddo! ITER = 1,NITER
1014 DO ii = 1,nindxglob
1015 i = indxglob(ii)
1016 pla(i) = pla(i)+ max(dpla1(i), zero)
1017 ENDDO
1018 ENDIF
1019 !========================================================================================
1020 ! LOCAL YIELD ALGORITHM
1021 !========================================================================================
1022 IF (nindxloc > 0) THEN
1023 DO ii = 1,nindxloc
1024 i = indxloc(ii)
1025 logz(i) = one
1026 logzm1(i) = one
1027 IF (uvar(i,3)/=zero)THEN
1028 psi2(i) = max(zero,(lnf2(i)*(teta2*pla1(i)))/(one+lnf2(i)) )
1029 phi2(i) = teta2*lnf2(i)/(one+lnf2(i))
1030 ENDIF
1031 IF (uvar(i,4)/=zero)THEN
1032 psi3(i) = max(zero,(lnf3(i)*(teta3*pla1(i)))/(one+lnf3(i)) )
1033 phi3(i) = teta3*lnf3(i)/(one+lnf3(i))
1034 ENDIF
1035 IF (uvar(i,5)/=zero)THEN
1036 psi4(i) = max(zero, (lnf4(i)*(teta4*pla1(i)))/(one+lnf4(i)) )
1037 phi4(i) = teta4*lnf4(i)/(one+lnf4(i))
1038 ENDIF
1039 IF (uvar(i,6)/=zero)THEN
1040 psi5(i) = max(zero,(lnf5(i)*(teta5*pla1(i)))/(one+lnf5(i)) )
1041 phi5(i) = teta5*lnf5(i)/(one+lnf5(i))
1042 ENDIF
1043 IF (totfrac(i) > zero)THEN
1044 logz(i) = log(totfrac(i))
1045 ENDIF
1046 IF ( totfracold(i)>zero)THEN !totfrac old
1047 logzm1(i) = log( totfracold(i))
1048 ENDIF
1049
1050 ! initialize incremental update
1051 stxx(i) = sxx(i)
1052 styy(i) = syy(i)
1053 stzz(i) = szz(i)
1054 stxy(i) = sxy(i)
1055 styz(i) = syz(i)
1056 stzx(i) = szx(i)
1057 svmt(i) = svm(i)
1058
1059
1060 a(i) = (one - totfrac(i))* gz(i) / e(i)
1061 b(i) = two * (alfa1 -alfa2)* (tempel(i)-tempo(i) ) *totfrac(i)*logz(i)
1062 c(i) = two*deth12(i)*abs((totfrac(i)*(one-logz(i))- totfracold(i)*(one-logzm1(i))))
1063
1064 hfct(i) = one + max(zero,seven_half*(svmo(i)/yld(i)-half) )
1065
1066 dpla(i) = a(i) * ( svm(i) - svmo(i) ) + b(i) + c(i) *hfct(i)
1067 gsig(i) = one/( one + three * (shear(i) / y1(i)) * dpla(i) )
1068
1069 normxx = three_half * soxx(i) /y1(i) ! DF/DSIG
1070 normyy = three_half * soyy(i) /y1(i)
1071 normzz = three_half * sozz(i) /y1(i)
1072 normxy = two * three_half * sigoxy(i) /y1(i)
1073 normyz = two * three_half * sigoyz(i) /y1(i)
1074 normzx = two * three_half * sigozx(i) /y1(i)
1075 dpxx(i) = dpla(i) * normxx
1076 dpyy(i) = dpla(i) * normyy
1077 dpzz(i) = dpla(i) * normzz
1078 dpxy(i) = dpla(i) * normxy
1079 dpyz(i) = dpla(i) * normyz
1080 dpzx(i) = dpla(i) * normzx
1081 ! Elasto-plastic stresses update
1082 signxx(i) = signxx(i) - (a1(i)*dpxx(i) + a2(i)*(dpyy(i) + dpzz(i)))
1083 signyy(i) = signyy(i) - (a1(i)*dpyy(i) + a2(i)*(dpxx(i) + dpzz(i)))
1084 signzz(i) = signzz(i) - (a1(i)*dpzz(i) + a2(i)*(dpxx(i) + dpyy(i)))
1085 signxy(i) = signxy(i) - shear(i)*dpxy(i)
1086 signyz(i) = signyz(i) - shear(i)*dpyz(i)
1087 signzx(i) = signzx(i) - shear(i)*dpzx(i)
1088 p(i) = -third*(signxx(i) + signyy(i) + signzz(i))
1089 ! Computation of the Von Mises equivalent stress of the stress tensor
1090 sxx(i) = signxx(i) + p(i)
1091 syy(i) = signyy(i) + p(i)
1092 szz(i) = signzz(i) + p(i)
1093 sxy(i) = signxy(i)
1094 syz(i) = signyz(i)
1095 szx(i) = signzx(i)
1096 svm(i) = half*(sxx(i)**2 + syy(i)**2 + szz(i)**2)
1097 . + (sxy(i)**2 + szx(i)**2 + syz(i)**2)
1098 svm(i) = sqrt(three*svm(i))
1099
1100 ENDDO !II = 1,NINDXLOC
1101 DO ii = 1,nindxloc
1102 i = indxloc(ii)
1103 dpla1(i)= dpla(i)/(one - totfrac(i))
1104 IF (uvar(i,3)>zero)THEN
1105 dpla2(i)= dpla1(i) *phi2(i) + psi2(i)
1106 ENDIF
1107 IF (uvar(i,4)>zero)THEN
1108 dpla3(i)= dpla1(i) *phi3(i) + psi3(i)
1109 ENDIF
1110 IF (uvar(i,5)>zero)THEN
1111 dpla4(i)= dpla1(i) *phi4(i) + psi4(i)
1112 ENDIF
1113 IF (uvar(i,6)>zero)THEN
1114 dpla5(i)= dpla1(i) *phi5(i) + psi5(i)
1115 ENDIF
1116 pla1(i)= uvar(i,17) + dpla1(i)
1117 pla2(i)= uvar(i,18) + dpla2(i)
1118 pla3(i)= uvar(i,19) + dpla3(i)
1119 pla4(i)= uvar(i,20) + dpla4(i)
1120 pla5(i)= uvar(i,21) + dpla5(i)
1121 !-------------------------------------
1122 !UPDATE yield for each phase
1123 !-------------------------------------
1124 xx1(i,1)=pla1(i)
1125 xx2(i,1)=pla2(i)
1126 xx3(i,1)=pla3(i)
1127 xx4(i,1)=pla4(i)
1128 xx5(i,1)=pla5(i)
1129 ipos1(i,1) = vartmp(i,1)
1130 ipos1(i,2) = vartmp(i,2)
1131 ipos1(i,3) = vartmp(i,3)
1132C
1133 ipos2(i,1) = vartmp(i,4)
1134 ipos2(i,2) = vartmp(i,5)
1135 ipos2(i,3) = vartmp(i,6)
1136
1137 ipos3(i,1) = vartmp(i,7)
1138 ipos3(i,2) = vartmp(i,8)
1139 ipos3(i,3) = vartmp(i,9)
1140C
1141 ipos4(i,1) = vartmp(i,10)
1142 ipos4(i,2) = vartmp(i,11)
1143 ipos4(i,3) = vartmp(i,12)
1144C
1145 ipos5(i,1) = vartmp(i,13)
1146 ipos5(i,2) = vartmp(i,14)
1147 ipos5(i,3) = vartmp(i,15)
1148 END DO
1149C
1150 CALL table_vinterp(table(itable(1)),nel,nel,ipos1,xx1,y1,dydx1)
1151 CALL table_vinterp(table(itable(2)),nel,nel,ipos2,xx2,y2,dydx2)
1152 CALL table_vinterp(table(itable(3)),nel,nel,ipos3,xx3,y3,dydx3)
1153 CALL table_vinterp(table(itable(4)),nel,nel,ipos4,xx4,y4,dydx4)
1154 CALL table_vinterp(table(itable(5)),nel,nel,ipos5,xx5,y5,dydx5)
1155 DO ii = 1,nindxloc
1156 i = indxloc(ii)
1157 y1(i)=y1(i)*yscale(1)
1158 y2(i)=y2(i)*yscale(2)
1159 y3(i)=y3(i)*yscale(3)
1160 y4(i)=y4(i)*yscale(4)
1161 y5(i)=y5(i)*yscale(5)
1162 dydx1(i)=dydx1(i)*yscale(1)
1163 dydx2(i)=dydx2(i)*yscale(2)
1164 dydx3(i)=dydx3(i)*yscale(3)
1165 dydx4(i)=dydx4(i)*yscale(4)
1166 dydx5(i)=dydx5(i)*yscale(5)
1167
1168 !new yield updated in newton iterations
1169 yld(i)=max(em20, y1(i)*frac1(i)+
1170 . y2(i)*frac2(i)+
1171 . y3(i)*frac3(i)+
1172 . y4(i)*frac4(i)+
1173 . y5(i)*frac5(i))
1174 yld(i)= yld(i)*faccs(i)
1175
1176 ENDDO !II = 1,NINDXLOC
1177 DO ii = 1,nindxloc
1178 i = indxloc(ii)
1179 vartmp(i,1) = ipos1(i,1)
1180 vartmp(i,2) = ipos1(i,2)
1181 vartmp(i,3) = ipos1(i,3)
1182C
1183 vartmp(i,4) = ipos2(i,1)
1184 vartmp(i,5) = ipos2(i,2)
1185 vartmp(i,6) = ipos2(i,3)
1186
1187 vartmp(i,7) = ipos3(i,1)
1188 vartmp(i,8) = ipos3(i,2)
1189 vartmp(i,9) = ipos3(i,3)
1190C
1191 vartmp(i,10) = ipos4(i,1)
1192 vartmp(i,11) = ipos4(i,2)
1193 vartmp(i,12) = ipos4(i,3)
1194C
1195 vartmp(i,13) = ipos5(i,1)
1196 vartmp(i,14) = ipos5(i,2)
1197 vartmp(i,15) = ipos5(i,3)
1198
1199 ENDDO !NINDXLOC
1200 ENDIF !(NINDXLOC > 0)
1201 DO i=1,nel
1202 IF (off(i) == one )THEN
1203 uvar(i,8) = tempel(i)
1204 uvar(i,17)= pla1(i) + max(dpla1(i), zero)
1205 uvar(i,18)= pla2(i) + max(dpla2(i), zero)
1206 uvar(i,19)= pla3(i) + max(dpla3(i), zero)
1207 uvar(i,20)= pla4(i) + max(dpla4(i), zero)
1208 uvar(i,21)= pla5(i) + max(dpla5(i), zero)
1209
1210
1211 uvar(i,35)= vol(i)
1212
1213 IF(uvar(i,15)<svm(i))uvar(i,15)=svm(i)
1214 IF(uvar(i,34)>tempel(i))uvar(i,34)=tempel(i)
1215
1216 uvar(i,23)= frac2(i)/xeq2 !X2(I)
1217 uvar(i,24)= frac3(i)/(one-xeq2) !X3(I)
1218 uvar(i,25)= frac4(i) !X4(I)
1219 uvar(i,44)= frac5(i)/max(em20,xgama(i))
1220
1221 uvar(i,2)= frac1(i)
1222 uvar(i,3)= frac2(i)
1223 uvar(i,4)= frac3(i)
1224 uvar(i,5)= frac4(i)
1225 uvar(i,6)= frac5(i)
1226 uvar(i,10)=xgama(i)
1227 hard(i)= zero
1228
1229 ENDIF ! off
1230 soundsp(i) = sqrt((rbulk(i)+four_over_3*shear(i))/rho0(i))
1231 uvar(i,9)=yld(i)
1232 ENDDO
1233 IF (alfa1/= zero .OR. alfa2 /= zero) THEN !compute internal thermal energy
1234 DO i=1,nel
1235 sigkk(i) = signxx(i)+signyy(i)+signzz(i)+sigoxx(i)+sigoyy(i)+sigozz(i)
1236 eintth(i) = eintth(i)-half*sigkk(i)*depsth(i)
1237 ENDDO
1238 ENDIF
1239 !IF(HEATFLAG == 0) THEN
1240 DO i=1,nel
1241 IF (off(i) == one )THEN
1242 IF (tempel(i)<ms.AND.tempel(i)<= uvar(i,8) .AND. tempel(i)<=1073.0)THEN !
1243 vr(i) = zero
1244 hard(i) = uvar(i,7)
1245 IF (uvar(i,16) >uvar(i,26))THEN
1246 vr(i) = (uvar(i,27)-uvar(i,33))*unitt/(uvar(i,16)-uvar(i,26))
1247 hard(i)= (frac2(i)+frac3(i))*hfp*log10(vr(i))+frac4(i)*hb+frac5(i)*hm
1248 uvar(i,7)= hard(i)
1249 ENDIF
1250 ENDIF
1251 die(i) =die(i)+ (lat2* (frac5(i)-uvar(i,6)) +
1252 . lat1* (frac2(i)-uvar(i,3)
1253 . + frac3(i)-uvar(i,4)
1254 . + frac4(i)-uvar(i,5) ) ) *vol(i)
1255
1256 ENDIF ! OFF(I) == ONE
1257 ENDDO
1258!-----------
1259 IF (jthe == 0) THEN
1260 temp(1:nel) = tempel(1:nel) ! for output
1261 END IF
1262!-----------
1263 RETURN
1264 END
alpha
#define alpha
Definition eval.h:35
interface_table_mod::table_vinterp
Definition table_mod.F:251
kirkaldykinetics
subroutine kirkaldykinetics(nel0, time, tempel, tempmin, ae1, ae3, bs, ms, fcfer, fcper, fcbai, fgrain, frac1, frac2, frac3, frac4, frac5, x2, x3, x4, x5, qr2, qr3, qr4, kper, kbain, alpha, xeq2, xeq4, xgama, totfrac, timestep, nicool, index, theaccfact)
Definition kirkaldykinetics.F:34
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
interface_table_mod
Definition table_mod.F:238
table_mod
Definition table_mod.F:112
phasekinetic2
subroutine phasekinetic2(nel0, time, tempel, tempo, tempmin, ae1, ae3, bs, ms, fcfer, fcper, fcbai, fgrain, frac1, frac2, frac3, frac4, frac5, x2, x3, x4, x5, gfac_f, phi_f, psi_f, cr_f, cf, gfac_p, phi_p, psi_p, cr_p, cp, gfac_b, phi_b, psi_b, cr_b, cb, phi_m, psi_m, n_m, fgfer, fgper, fgbai, qr2, qr3, qr4, kper, kbain, alpha, xeq2, xeq4, xgama, totfrac, timestep, nicool, index, theaccfact)
Definition phasekinetic2.F:35
sigeps80
subroutine sigeps80(nel, nuparam, nuvar, mfunc, kfunc, npf, tf, time, timestep, uparam, rho0, rho, vol, eint, 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, soundsp, viscmax, uvar, off, ngl, ipm, mat, epsp, yld, pla, table, temp, nvartmp, vartmp, trdepsth, eintth, jthe, idt_therm, theaccfact)
Definition sigeps80.F:53
table_mod::ttable
Definition table_mod.F:125