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