OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps71c.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/.
23C New law -- material shape memory
24!||====================================================================
25!|| sigeps71c ../engine/source/materials/mat/mat071/sigeps71c.F
26!||--- called by ------------------------------------------------------
27!|| mulawc ../engine/source/materials/mat_share/mulawc.F90
28!||--- calls -----------------------------------------------------
29!|| finter ../engine/source/tools/curve/finter.F
30!||====================================================================
31 SUBROUTINE sigeps71c(
32 1 NEL, NUPARAM, NUVAR, MFUNC,
33 2 KFUNC, NPF, NPT0, IPT,
34 3 IFLAG, TF, TIME, TIMESTEP,
35 4 UPARAM, RHO0, AREA, EINT,
36 5 THKLY, EPSPXX, EPSPYY, EPSPXY,
37 6 EPSPYZ, EPSPZX, DEPSXX, DEPSYY,
38 7 DEPSXY, DEPSYZ, DEPSZX, EPSXX,
39 8 EPSYY, EPSXY, EPSYZ, EPSZX,
40 9 SIGOXX, SIGOYY, SIGOXY, SIGOYZ,
41 A SIGOZX, SIGNXX, SIGNYY, SIGNXY,
42 B SIGNYZ, SIGNZX, SIGVXX, SIGVYY,
43 C SIGVXY, SIGVYZ, SIGVZX, SOUNDSP,
44 D VISCMAX, THK, PLA, UVAR,
45 E OFF, NGL, IPM, MAT,
46 F ETSE, GS, YLD, VOL,
47 G TEMP, ISMSTR, JTHE)
48C-----------------------------------------------
49C I m p l i c i t T y p e s
50C-----------------------------------------------
51#include "implicit_f.inc"
52C-----------------------------------------------
53C G l o b a l P a r a m e t e r s
54C-----------------------------------------------
55#include "mvsiz_p.inc"
56C---------+---------+---+---+--------------------------------------------
57C VAR | SIZE |TYD| RW| DEFINITION
58C---------+---------+---+---+--------------------------------------------
59C NEL0 | 1 | I | R | SIZE OF THE ELEMENT GROUP NEL0
60C NUPARAM | 1 | I | R | SIZE OF THE USER PARAMETER ARRAY
61C NUVAR | 1 | I | R | NUMBER OF USER ELEMENT VARIABLES
62C---------+---------+---+---+--------------------------------------------
63C NFUNC | 1 | I | R | NUMBER FUNCTION USED FOR THIS USER LAW
64C IFUNC | NFUNC | I | R | FUNCTION INDEX
65C NPF | * | I | R | FUNCTION ARRAY
66C NPT0 | 1 | I | R | NUMBER OF LAYERS OR INTEGRATION POINTS
67C IPT | 1 | I | R | LAYER OR INTEGRATION POINT NUMBER
68C IFLAG | * | I | R | GEOMETRICAL FLAGS
69C TF | * | F | R | FUNCTION ARRAY
70C---------+---------+---+---+--------------------------------------------
71C TIME | 1 | F | R | CURRENT TIME
72C TIMESTEP| 1 | F | R | CURRENT TIME STEP
73C UPARAM | NUPARAM | F | R | USER MATERIAL PARAMETER ARRAY
74C RHO0 | NEL0 | F | R | INITIAL DENSITY
75C AREA | NEL0 | F | R | AREA
76C EINT | 2*NEL0 | F | R | INTERNAL ENERGY(MEMBRANE,BENDING)
77C THKLY | NEL0 | F | R | LAYER THICKNESS
78C EPSPXX | NEL0 | F | R | STRAIN RATE XX
79C EPSPYY | NEL0 | F | R | STRAIN RATE YY
80C ... | | | |
81C DEPSXX | NEL0 | F | R | STRAIN INCREMENT XX
82C DEPSYY | NEL0 | F | R | STRAIN INCREMENT YY
83C ... | | | |
84C EPSXX | NEL0 | F | R | STRAIN XX
85C EPSYY | NEL0 | F | R | STRAIN YY
86C ... | | | |
87C SIGOXX | NEL0 | F | R | OLD ELASTO PLASTIC STRESS XX
88C SIGOYY | NEL0 | F | R | OLD ELASTO PLASTIC STRESS YY
89C ... | | | |
90C---------+---------+---+---+--------------------------------------------
91C SIGNXX | NEL0 | F | W | NEW ELASTO PLASTIC STRESS XX
92C SIGNYY | NEL0 | F | W | NEW ELASTO PLASTIC STRESS YY
93C ... | | | |
94C SIGVXX | NEL0 | F | W | VISCOUS STRESS XX
95C SIGVYY | NEL0 | F | W | VISCOUS STRESS YY
96C ... | | | |
97C SOUNDSP | NEL0 | F | W | SOUND SPEED (NEEDED FOR TIME STEP)
98C VISCMAX | NEL0 | F | W | MAXIMUM DAMPING MODULUS(NEEDED FOR TIME STEP)
99C---------+---------+---+---+--------------------------------------------
100C THK | NEL0 | F |R/W| THICKNESS
101C PLA | NEL0 | F |R/W| PLASTIC STRAIN
102C UVAR |NEL0*NUVAR| F |R/W| USER ELEMENT VARIABLE ARRAY
103C OFF | NEL0 | F |R/W| DELETED ELEMENT FLAG (=1. ON, =0. OFF)
104C---------+---------+---+---+--------------------------------------------
105C C o m m o n B l o c k s
106C-----------------------------------------------
107#include "param_c.inc"
108C-----------------------------------------------
109C I N P U T A r g u m e n t s
110C-----------------------------------------------
111C
112 INTEGER, INTENT(IN) :: ISMSTR, JTHE
113 INTEGER NEL, NUPARAM, NUVAR, NPT0, IPT,IFLAG(*),
114 . NGL(NEL),MAT(NEL),IPM(NPROPMI,*)
115 my_real TIME,TIMESTEP,UPARAM(*),
116 . AREA(NEL),RHO0(NEL),EINT(NEL,2),
117 . THKLY(NEL),PLA(NEL),
118 . EPSPXX(NEL),EPSPYY(NEL),
119 . EPSPXY(NEL),EPSPYZ(NEL),EPSPZX(NEL),
120 . DEPSXX(NEL),DEPSYY(NEL),
121 . DEPSXY(NEL),DEPSYZ(NEL),DEPSZX(NEL),
122 . EPSXX(NEL) ,EPSYY(NEL) ,
123 . EPSXY(NEL) ,EPSYZ(NEL) ,EPSZX(NEL) ,
124 . SIGOXX(NEL),SIGOYY(NEL),
125 . SIGOXY(NEL),SIGOYZ(NEL),SIGOZX(NEL),
126 . GS(*) ,VOL(NEL) , TEMP(NEL)
127C-----------------------------------------------
128C O U T P U T A r g u m e n t s
129C-----------------------------------------------
130 my_real
131 . signxx(nel),signyy(nel),
132 . signxy(nel),signyz(nel),signzx(nel),
133 . sigvxx(nel),sigvyy(nel),
134 . sigvxy(nel),sigvyz(nel),sigvzx(nel),
135 . soundsp(nel),viscmax(nel),etse(nel)
136C-----------------------------------------------
137C I N P U T O U T P U T A r g u m e n t s
138C-----------------------------------------------
139 my_real uvar(nel,nuvar), off(nel),thk(nel),yld(nel)
140C-----------------------------------------------
141C VARIABLES FOR FUNCTION INTERPOLATION
142C-----------------------------------------------
143 INTEGER NPF(*), MFUNC, KFUNC(MFUNC)
144 my_real FINTER ,TF(*)
145 EXTERNAL FINTER
146C Y = FINTER(IFUNC(J),X,NPF,TF,DYDX)
147C Y : y = f(x)
148C X : x
149C DYDX : f'(x) = dy/dx
150C IFUNC(J): FUNCTION INDEX
151C J : FIRST(J=1), SECOND(J=2) .. FUNCTION USED FOR THIS LAW
152C NPF,TF : FUNCTION PARAMETER
153C-----------------------------------------------
154C L o c a l V a r i a b l e s
155C-----------------------------------------------
156 INTEGER I,JJ,KK,ITERK,EFLAG
157 INTEGER CRIT_LOOP(NEL)
158 INTEGER INDX_LOOP(NEL),NINDX_LOOP,NINDX_LOOP_LOC,II
159 INTEGER, DIMENSION(NEL) :: INDX_LOOP_LOC
160 .
161 my_real
162 . e,emart,nu,g,g2,sqdt,a,b,c,fct,fctp, dftr,unmxn,db,
163 . alpha,epsl,aa1,fm,dfmsa,dfmas,beta,gm,km,h,
164 . dgt,dkt,cp,tini,sspsh ,sspsol,
165 . k,p11(nel),sxx(nel),cas,csa,tsas,tfas, tssa,tfsa,inve,
166 . syy(nel),szz(nel),fass,fsas,fasf,fsaf,rsas,rfas,
167 . sv(nel),fs(nel),fs0,yld_ass,yld_asf,yld_sas,yld_saf,rssa,rfsa,
168 . dfm(nel),rv_pui,
169 . test,
170 . ne,
171 . nxx(nel),nyy(nel),nzz(nel),
172 . trde(nel),
173 . gt(nel),kt(nel),
174 . ee1(nel),ee2(nel),ee3(nel),det(nel),
175 . sigxx(nel), sigyy(nel), sigzz(nel),dfs(nel)
176 my_real
177 . epszz(nel),
178 . sigzim1(nel),epsim1(nel),epsi(nel)
179 my_real
180 . ev(mvsiz,3)
181 my_real
182 . eigv(nel,3,2),trav(nel),rootv(nel),evv(nel,3)
183
184 DATA iterk/10/
185C======================================================================|
186 e = uparam(1)
187 nu = uparam(2)
188 g = uparam(3)
189 k = uparam(4)
190 aa1 = uparam(5)
191 yld_ass = uparam(6)
192 yld_asf = uparam(7)
193 yld_sas = uparam(8)
194 yld_saf = uparam(9)
195 alpha = uparam(10)
196 epsl = uparam(11)/(sqrt(two_third)+alpha)
197 emart = uparam(14)
198 eflag = int(uparam(15))
199 gm = uparam(16)
200 km = uparam(17)
201 g2 = two*g
202 beta = epsl*(g2+nine*k*alpha*alpha)
203 sqdt = sqrt(two/three)
204 cas = uparam(18)
205 csa = uparam(19)
206 tsas = uparam(20)
207 tfas = uparam(21)
208 tssa = uparam(22)
209 tfsa = uparam(23)
210 cp = uparam(24)
211 tini = uparam(25)
212
213C principal stretch (def gradient eigenvalues)
214 trav(1:nel) = epsxx(1:nel)+epsyy(1:nel)
215 rootv(1:nel) = sqrt((epsxx(1:nel)-epsyy(1:nel))*(epsxx(1:nel)-epsyy(1:nel))
216 . + epsxy(1:nel)*epsxy(1:nel))
217 evv(1:nel,1) = half*(trav(1:nel)+rootv(1:nel))
218 evv(1:nel,2) = half*(trav(1:nel)-rootv(1:nel))
219 evv(1:nel,3) = zero
220C rot matrix (eigenvectors)
221 DO i=1,nel
222 IF(abs(evv(i,2)-evv(i,1)) < em10) THEN
223 eigv(i,1,1) = one
224 eigv(i,2,1) = one
225 eigv(i,3,1) = zero
226 eigv(i,1,2) = zero
227 eigv(i,2,2) = zero
228 eigv(i,3,2) = zero
229 ELSE
230 eigv(i,1,1) = (epsxx(i)-evv(i,2)) /rootv(i)
231 eigv(i,2,1) = (epsyy(i)-evv(i,2)) /rootv(i)
232 eigv(i,1,2) = (evv(i,1)-epsxx(i)) /rootv(i)
233 eigv(i,2,2) = (evv(i,1)-epsyy(i)) /rootv(i)
234 eigv(i,3,1) = (half*epsxy(i)) /rootv(i)
235 eigv(i,3,2) =-(half*epsxy(i)) /rootv(i)
236 ENDIF
237 ENDDO
238 IF (ismstr == 1 .OR. ismstr == 3 .OR. ismstr == 11) THEN ! engineering strain
239 DO i=1,nel
240 ev(i,1)=evv(i,1)+ one
241 ev(i,2)=evv(i,2)+ one
242 ev(i,3)=one/ev(i,1)/ev(i,2)
243 ENDDO
244 ELSEIF(ismstr == 10) THEN
245 DO i=1,nel
246 ev(i,1)=sqrt(evv(i,1)+ one)
247 ev(i,2)=sqrt(evv(i,2)+ one)
248 ev(i,3)=one/ev(i,1)/ev(i,2)!initialization before plane stress iterations
249 ENDDO
250 ELSE ! true strain
251 DO i=1,nel
252 ev(i,1)=exp(evv(i,1))
253 ev(i,2)=exp(evv(i,2))
254 ev(i,3)=one/ev(i,1)/ev(i,2)
255 ENDDO
256 ENDIF
257C-----------------------------------------------
258C calculation of the temperature (conduction or adiabatic)
259C--------------------
260 IF (jthe == 0) THEN
261 DO i=1,nel
262 temp(i) = tini + (eint(i,1)+ eint(i,2))/vol(i)/ rho0(i)/cp
263 ENDDO
264 ENDIF
265 !compute limits for start and end of transformation
266 rsas = yld_ass *(sqdt+alpha)-cas*tsas !start from aust to martensite
267 rfas = yld_asf *(sqdt+alpha)-cas*tfas !end from aust to martensite
268 rssa = yld_sas *(sqdt+alpha)-csa*tssa !start from martensite to austenite
269 rfsa = yld_saf *(sqdt+alpha)-csa*tfsa !end from martensite to austenite
270 !----------------------------------------------
271 crit_loop(1:nel) = 0
272 DO i=1,nel
273 indx_loop(i) = i
274 indx_loop_loc(i) = 0
275 ENDDO
276 nindx_loop = nel
277c----------------------------------------------
278c----------------------------------------------
279c DEBUT BOUCLE 1 SUR K plane stress
280c----------------------------------------------
281c----------------------------------------------
282 DO kk=1,iterk ! plane stress iteration
283#include "vectorize.inc"
284 DO ii =1,nindx_loop
285 i = indx_loop(ii)
286 !DET A (A = GRADIENT OF DEFORMATION)
287 det(i) =ev(i,1)*ev(i,2)*ev(i,3)
288 IF(det(i)/=zero) THEN
289 trde(i) = log(det(i)) !TRACE OF DEFORMATION
290 rv_pui = exp((-third)*trde(i))
291 ELSE
292 rv_pui = zero
293 trde(i)= zero
294 ENDIF
295 ee1(i) = log(ev(i,1)*rv_pui) !DEVIATOR OF DEFORMATION
296 ee2(i) = log(ev(i,2)*rv_pui)
297 ee3(i) = log(ev(i,3)*rv_pui)
298 ENDDO
299 !==============================================================
300 IF (eflag > zero)THEN !dependency on martensite fraction
301 !==============================================================
302#include "vectorize.inc"
303 DO ii =1,nindx_loop
304 i = indx_loop(ii)
305 fm = uvar(i,1) ! fraction of Martensite
306 gt(i) = g + fm * (gm - g) !G_n
307 kt(i) = k + fm * (km - k) !K_n
308 !pressure
309 p11(i) = kt(i) * (trde(i) - three*alpha*epsl*fm)
310 ! n= e/ norm(e)
311 ne = sqrt( ee1(i)**2 + ee2(i)**2 + ee3(i)**2 )
312 nxx(i) =ee1(i)/max(ne,em20)
313 nyy(i) =ee2(i)/max(ne,em20)
314 nzz(i) =ee3(i)/max(ne,em20)
315
316! Estimation dev(sigma_n+1)
317 sxx(i)= two*gt(i)*(ee1(i) -epsl*fm*nxx(i))
318 syy(i)= two*gt(i)*(ee2(i) -epsl*fm*nyy(i))
319 szz(i)= two*gt(i)*(ee3(i) -epsl*fm*nzz(i))
320
321 sv(i) = sqrt( sxx(i)*sxx(i) + syy(i)*syy(i) + szz(i)*szz(i) )
322 !-------------------
323 ! transformation
324 !-------------------
325 dfmsa = zero
326 dfmas = zero
327
328 ! loading function
329 fs(i) = sv(i) + three*alpha*p11(i) - cas*temp(i) ! loading function A--> M
330 !----------------------
331 ! Check Austenite -----> martensite
332 fass = fs(i) -rsas
333 fasf = fs(i) -rfas
334 fs0 = uvar(i,2)
335 beta = epsl*(two*gt(i)+nine*kt(i)*alpha*alpha)
336 IF((fs(i) - fs0) > zero .AND. fass > zero.AND. fasf < zero .AND. fm < one )THEN
337 ! DFMAS = MIN(ONE, -(FS-FS0)*(ONE-FM)/(FASF-BETA*(ONE-FM) ) ) ! (should be positive)
338 db = (two * (gm-g) +nine*alpha*alpha*(km-k)) *epsl
339 unmxn = one - fm
340 dftr = two*ne*(gm-g) + three*alpha*trde(i)*(km-k)
341 dfmas = min(one, -(fs(i)-fs0)*(one-fm)/(fasf-beta*(one-fm) ) )
342 a = unmxn *db
343 b = (rfas-fs(i)+unmxn*(beta-dftr))
344 c = unmxn*(fs0 - fs(i))
345 DO jj = 1,3
346 fct = dfmas*dfmas *a+ dfmas* b +c
347 fctp = two*dfmas *a+ b
348 dfmas = dfmas - fct / fctp
349 ENDDO
350 dfmas = min(one,dfmas ) ! (should be positive)
351 ENDIF
352
353
354 !IF((FS(I) - FS0) > ZERO .AND. FASS > ZERO.AND. FASF < ZERO.AND. FM < ONE ) THEN!
355 ! DFMAS = MIN(ONE, -(FS(I)-FS0)*(ONE-FM)/(FASF-BETA*(ONE-FM) ) ) ! (should be positive)
356 !ENDIF
357 !----------------------
358 ! Check marteniste -----> austenite (dependance au tempurature a voir ...)
359 !----------------------
360 fs(i) = sv(i) + three*alpha*p11(i) - csa*temp(i) ! Unloading function M--> A
361 fsas = fs(i) - rssa
362 fsaf = fs(i) - rfsa
363 fs0 = uvar(i,3)
364 IF((fs(i) - fs0) < zero .AND. fsas < zero.AND. fsaf > zero )THEN
365 !DFMSA = MAX(-ONE , FM * (FS - FS0)/ (FSAF+BETA*FM) )
366 db = (two * (gm-g) +nine*alpha*alpha*(km-k))*epsl
367 dftr = two * (gm-g)*ne+ three*alpha*(km-k)*trde(i)
368 dfmsa = zero
369 a = fm *db
370 b = -(rfsa-fs(i)+fm*(dftr-beta))
371 c = -fm*(fs(i) - fs0)
372 DO jj = 1,3
373 fct = dfmsa*dfmsa *a+ dfmsa* b +c
374 fctp = two*dfmsa *a+ b
375 dfmsa = dfmsa - fct / fctp
376 ENDDO
377 dfmsa = max(-one , dfmsa )
378 ENDIF
379 !IF((FS(I) - FS0) < ZERO .AND. FSAS < ZERO .AND. FSAF > ZERO ) THEN
380 ! DFMSA = MAX(-ONE , FM * (FS(I) - FS0)/ (FSAF+BETA*FM) ) !compute marteniste fraction
381 !ENDIF
382
383
384 dfm(i) = dfmas + dfmsa ! new martensite fraction
385 IF(dfm(i) < zero .AND. fm == zero) dfm(i) = zero
386c
387
388 !UPDATE
389 dgt = dfm(i) * (gm - g)
390 dkt = dfm(i) * (km - k)
391 sxx(i) = sxx(i) -two*gt(i)* epsl*dfm(i)*nxx(i)
392 . + two*dgt* (ee1(i)-epsl*nxx(i)*dfm(i))
393 syy(i) = syy(i) -two*gt(i)* epsl*dfm(i)*nyy(i)
394 . + two*dgt* (ee2(i)-epsl*nyy(i)*dfm(i))
395 szz(i) = szz(i) -two*gt(i)* epsl*dfm(i)*nzz(i)
396 . + two*dgt* (ee3(i)-epsl*nzz(i)*dfm(i))
397
398 p11(i) = p11(i) - kt(i)*epsl*three*alpha*dfm(i)
399 . + dkt *(trde(i) -epsl*three*alpha*dfm(i))
400
401 !KIRCHHOFF STRESS
402 sigxx(i)= sxx(i) + p11(i)
403 sigyy(i)= syy(i) + p11(i)
404 sigzz(i)= szz(i) + p11(i)
405
406 sv(i) = sqrt( sxx(i)*sxx(i) + syy(i)*syy(i) + szz(i)*szz(i) )
407 fs(i) = sv(i) + three*alpha*p11(i)
408 dfs(i) = zero
409 IF (dfmas /= zero) dfs(i) = abs(fs(i)-cas*temp(i) - fs0)
410 IF (dfmsa /= zero) dfs(i) = abs(fs(i)-csa*temp(i) - fs0)
411 IF (dfs(i) /= zero .AND. epsl /= zero .AND. dfm(i)/= zero) THEN
412 h = dfs(i)/epsl/dfm(i)
413 !ETSE(I) = H / (H+YOUNG) for shells
414 etse(i) = h /( (e + uvar(i,1)*(emart-e)) + h)
415 ELSE
416 etse(i) = one
417 ENDIF
418 ENDDO
419 ! ----------------
420#include "vectorize.inc"
421 DO ii =1,nindx_loop
422 i = indx_loop(ii)
423 IF(det(i)/=zero)THEN
424 inve = one/det(i)
425 ELSE
426 inve = zero
427 ENDIF
428 !CAUCHY STRESS
429 sigxx(i)= sigxx(i) *inve
430 sigyy(i)= sigyy(i) *inve
431 sigzz(i)= sigzz(i) *inve
432 ENDDO
433 !==============================================================
434 ELSE ! no dependency martensite fraction
435 !==============================================================
436#include "vectorize.inc"
437 DO ii =1,nindx_loop
438 i = indx_loop(ii)
439 fm = uvar(i,1) ! fraction of Martensite
440 !pressure
441 p11(i) = k * (trde(i) - three*alpha*epsl*fm)
442 ! n= e/ norm(e)
443 ne = sqrt( ee1(i)**2 + ee2(i)**2 + ee3(i)**2 )
444 nxx(i) =ee1(i)/(ne+em20)
445 nyy(i) =ee2(i)/(ne+em20)
446 nzz(i) =ee3(i)/(ne+em20)
447
448 !Estimation dev(sigma_n+1)
449 sxx(i)= g2*(ee1(i) -epsl*fm*nxx(i))
450 syy(i)= g2*(ee2(i) -epsl*fm*nyy(i))
451 szz(i)= g2*(ee3(i) -epsl*fm*nzz(i))
452 sv(i) = sqrt( sxx(i)*sxx(i) + syy(i)*syy(i) + szz(i)*szz(i) )
453
454 !-------------------
455 ! transformation
456 !-------------------
457 dfmsa = zero
458 dfmas = zero
459
460 ! loading function
461 fs(i) = sv(i) + three*alpha*p11(i) - cas*temp(i) ! loading function A--> M
462 !----------------------
463 !Check Austenite -----> Martensite (Depending Temperature to see ...)
464 fass = fs(i) -rsas
465 fasf = fs(i) -rfas
466 fs0 = uvar(i,2)
467 IF((fs(i) - fs0) > zero .AND. fass > zero.AND. fasf < zero.AND. fm < one ) then!
468 dfmas = min(one, -(fs(i)-fs0)*(one-fm)/(fasf-beta*(one-fm) ) ) ! (should be positive)
469 ENDIF
470 !----------------------
471 !Check marteniste -----> austenite (dependance au tempurature a voir ...)
472 fs(i) = sv(i) + three*alpha*p11(i) - csa*temp(i) ! Unloading function M--> A
473 fsas = fs(i) - rssa
474 fsaf = fs(i) - rfsa
475 fs0 = uvar(i,3)
476 IF((fs(i) - fs0) < zero .AND. fsas < zero .AND. fsaf > zero ) THEN
477 dfmsa = max(-one , fm * (fs(i) - fs0)/ (fsaf+beta*fm) ) !compute marteniste fraction
478 ENDIF
479
480 dfm(i) = dfmas + dfmsa !new martensite fraction
481 IF(dfm(i) < zero .AND. fm == zero) dfm(i) = zero
482
483 !UPDATE
484 sxx(i) = sxx(i) -g2* epsl*dfm(i)*nxx(i)
485 syy(i) = syy(i) -g2* epsl*dfm(i)*nyy(i)
486 szz(i) = szz(i) -g2* epsl*dfm(i)*nzz(i)
487 p11(i) = p11(i) - k*epsl*three*alpha*dfm(i)
488 sigxx(i)= sxx(i) + p11(i)
489 sigyy(i)= syy(i) + p11(i)
490 sigzz(i)= szz(i) + p11(i)
491 sv(i) = sqrt( sxx(i)*sxx(i) + syy(i)*syy(i) + szz(i)*szz(i) )
492 fs(i) = sv(i) + three*alpha*p11(i)
493 dfs(i) = zero
494 IF (dfmas /= zero) dfs(i) = abs(fs(i)-cas*temp(i) - fs0)
495 IF (dfmsa /= zero) dfs(i) = abs(fs(i)-csa*temp(i) - fs0)
496 IF (dfs(i) /= zero .AND. epsl /= zero.AND.dfm(i)/= zero) THEN
497 h = dfs(i)/epsl/dfm(i)
498 etse(i) = h /(e + h)
499 ELSE
500 etse(i) = one
501 ENDIF
502 ENDDO
503 ! ----------------
504#include "vectorize.inc"
505 DO ii =1,nindx_loop
506 i = indx_loop(ii)
507 IF(det(i)/=zero)THEN
508 inve = one/det(i)
509 ELSE
510 inve = zero
511 ENDIF
512 !CAUCHY STRESS
513 sigxx(i)= sigxx(i) *inve
514 sigyy(i)= sigyy(i) *inve
515 sigzz(i)= sigzz(i) *inve
516 ENDDO
517
518 !==============================================================
519 ENDIF ! dependency on martensite fraction
520 !==============================================================
521
522 !-----------------------------------------------------
523 !check plane stress condition
524 !-----------------------------------------------------
525 nindx_loop_loc = 0
526 DO ii =1,nindx_loop
527 i = indx_loop(ii)
528 IF(abs(sigzz(i))>em20.OR.kk< 3)THEN ! CRITERE
529 IF (kk == 1) THEN
530 epsim1(i) = ev(i,3)
531 ev(i,3) = ev(i,3) /two
532 sigzim1(i) = sigzz(i)
533 ELSE
534 test = sigzz(i)-sigzim1(i)
535 epsi(i) = ev(i,3)
536 IF (test/=zero)THEN
537 ev(i,3) = ev(i,3)-sigzz(i) *(ev(i,3)-epsim1(i))/
538 . (sigzz(i)-sigzim1(i))
539 ELSE
540 ev(i,3) = ev(i,3)-sigzz(i) *(ev(i,3)-epsim1(i))/
541 . em10
542 ENDIF
543 epsim1(i) = epsi(i)
544 sigzim1(i) = sigzz(i)
545 !COMPUTE SECANT EXPRESSION
546 ENDIF
547 nindx_loop_loc = nindx_loop_loc + 1
548 indx_loop_loc(nindx_loop_loc) = i
549 ENDIF
550 ENDDO ! 1:NEL
551 ! ----------------
552 nindx_loop = nindx_loop_loc
553 indx_loop(1:nindx_loop) = indx_loop_loc(1:nindx_loop_loc)
554
555 ! -------------- -------------------------------
556 ENDDO ! plane stress iteration
557 ! ----------------------------------------------
558
559
560C ---------------update UVAR -------------------------------
561 uvar(1:nel,1) = uvar(1:nel,1) + dfm(1:nel)
562 uvar(1:nel,1) = max(zero, uvar(1:nel,1))
563 uvar(1:nel,2) = fs(1:nel)- cas*temp(1:nel)
564 uvar(1:nel,3) = fs(1:nel)- csa*temp(1:nel)
565C ----------------------------------------------------------
566 DO i=1,nel
567 epszz(i) =ev(i,3) - one ! left gauchy-green strain
568 ENDDO
569C TRANSFORM PRINCIPAL CAUCHY STRESSES TO GLOBAL DIRECTIONS
570 signxx(1:nel) = eigv(1:nel,1,1)*sigxx(1:nel)+eigv(1:nel,1,2)*sigyy(1:nel)
571 signyy(1:nel) = eigv(1:nel,2,1)*sigxx(1:nel)+eigv(1:nel,2,2)*sigyy(1:nel)
572 signxy(1:nel) = eigv(1:nel,3,1)*sigxx(1:nel)+eigv(1:nel,3,2)*sigyy(1:nel)
573 signyz(1:nel) = sigoyz(1:nel) + gs(1:nel)*depsyz(1:nel)
574 signzx(1:nel) = sigozx(1:nel) + gs(1:nel)*depszx(1:nel)
575 thk(1:nel) = thk(1:nel) + (epszz(1:nel)-uvar(1:nel,4))*thkly(1:nel)*off(1:nel)
576
577C ---------------update UVAR ---------------------
578 uvar(1:nel,1) = min(one, uvar(1:nel,1))
579 uvar(1:nel,10) = epsxx(1:nel)
580 uvar(1:nel,4) = epszz(1:nel)
581 uvar(1:nel,8) = sigzz(1:nel)
582c sound velocity
583 viscmax(1:nel) = zero
584 DO i=1,nel
585 sspsh = sqrt( e /(one - nu*nu)/rho0(i) )
586 sspsol = sqrt( aa1/rho0(i) )
587 soundsp(i) =max(sspsh ,sspsol)
588 ENDDO
589C------------------------------------------
590 RETURN
591 END SUBROUTINE sigeps71c
alpha
#define alpha
Definition eval.h:35
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
sigeps71c
subroutine sigeps71c(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, ipm, mat, etse, gs, yld, vol, temp, ismstr, jthe)
Definition sigeps71c.F:48