OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps73c.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!|| sigeps73c ../engine/source/materials/mat/mat073/sigeps73c.F
25!||--- called by ------------------------------------------------------
26!|| mulawc ../engine/source/materials/mat_share/mulawc.F90
27!||--- calls -----------------------------------------------------
28!|| finter ../engine/source/tools/curve/finter.F
29!|| table_vinterp ../engine/source/tools/curve/table_tools.F
30!||--- uses -----------------------------------------------------
31!|| interface_table_mod ../engine/share/modules/table_mod.F
32!|| table_mod ../engine/share/modules/table_mod.F
33!||====================================================================
34 SUBROUTINE sigeps73c(
35 1 NEL, NUPARAM, NUVAR, TIME,
36 2 TIMESTEP, UPARAM, RHO0, AREA,
37 3 EINT, THKLY, EPSPXX, EPSPYY,
38 4 EPSPXY, EPSPYZ, EPSPZX, DEPSXX,
39 5 DEPSYY, DEPSXY, DEPSYZ, DEPSZX,
40 6 EPSXX, EPSYY, EPSXY, EPSYZ,
41 7 EPSZX, SIGOXX, SIGOYY, SIGOXY,
42 8 SIGOYZ, SIGOZX, SIGNXX, SIGNYY,
43 9 SIGNXY, SIGNYZ, SIGNZX, SIGVXX,
44 A SIGVYY, SIGVXY, SIGVYZ, SIGVZX,
45 B SOUNDSP, VISCMAX, THK, PLA,
46 C UVAR, OFF, NGL, ITABLE,
47 D ETSE, GS, SIGY, DPLA1,
48 E EPSP, TABLE, VOL, TEMPEL,
49 F DIE, NPF, NFUNC,
50 G IFUNC, TF, SHF, HARDM,
51 H SEQ_OUTPUT,INLOC, DPLANL, JTHE,
52 I LOFF)
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 G L O B A L P A R A M E T E R S
62C-----------------------------------------------
63#include "mvsiz_p.inc"
64C-----------------------------------------------
65C I N P U T A R G U M E N T S
66C-----------------------------------------------
67 INTEGER, INTENT(IN) :: JTHE
68 INTEGER NEL, NUPARAM, NUVAR, NGL(NEL),ITABLE(*),
69 . INLOC
70 my_real
71 . TIME,TIMESTEP,UPARAM(NUPARAM),
72 . area(nel),rho0(nel),eint(nel,2),
73 . thkly(nel),pla(nel),
74 . epspxx(nel),epspyy(nel),
75 . epspxy(nel),epspyz(nel),epspzx(nel),
76 . depsxx(nel),depsyy(nel),
77 . depsxy(nel),depsyz(nel),depszx(nel),
78 . epsxx(nel) ,epsyy(nel) ,
79 . epsxy(nel) ,epsyz(nel) ,epszx(nel) ,
80 . sigoxx(nel),sigoyy(nel),
81 . sigoxy(nel),sigoyz(nel),sigozx(nel),shf(nel),
82 . gs(nel),vol(nel),tempel(nel),die(nel),hardm(*),
83 . seq_output(nel),dplanl(nel)
84 TYPE(ttable) TABLE(*)
85 my_real, DIMENSION(NEL), INTENT(IN) :: loff
86C-----------------------------------------------
87C O U T P U T A R G U M E N T S
88C-----------------------------------------------
89 my_real
90 . signxx(nel),signyy(nel),
91 . signxy(nel),signyz(nel),signzx(nel),
92 . sigvxx(nel),sigvyy(nel),sigy(nel),
93 . sigvxy(nel),sigvyz(nel),sigvzx(nel),
94 . soundsp(nel),viscmax(nel),etse(nel)
95C-----------------------------------------------
96C I N P U T O U T P U T A R G U M E N T S
97C-----------------------------------------------
98 my_real uvar(nel,nuvar), off(nel),thk(nel),
99 . dpla1(nel),epsp(nel)
100C-----------------------------------------------
101C VARIABLES FOR FUNCTION INTERPOLATION
102C-----------------------------------------------
103 INTEGER NPF(*), NFUNC, IFUNC(NFUNC)
104 my_real
105 . FINTER ,TF(*)
106 EXTERNAL FINTER
107C-----------------------------------------------
108C L O C A L V A R I A B L E S
109C-----------------------------------------------
110 INTEGER I,J,N,NINDX,NMAX,
111 . INDEX(MVSIZ),OPTE,IFUNCE,FUNC_TAB
112 my_real
113 . A,B,S12,DEZZ,S1,S2,S3,EPST,
114 . A_1,A_2,A_3,P_1,P_2,P_3,PP1,PP2,PP3,DR,DR0,
115 . E,A1,A2,G,G3,
116 . a01,a02,a03,a12,
117 . axx(mvsiz),ayy(mvsiz),axy(mvsiz),a_xy(mvsiz),
118 . dydx(mvsiz),svm(mvsiz),
119 . yld(mvsiz),yk(mvsiz),
120 . nnu1,nu,nu2,nu3,nu4,nu5,
121 . jq(mvsiz),jq2(mvsiz),
122 . epsmax,dpla_i(mvsiz),dpla_j(mvsiz),fail(mvsiz),
123 . epsr1,epsr2,s11(mvsiz),s22(mvsiz),
124 . b_1(mvsiz),b_2(mvsiz),b_3(mvsiz),q12(mvsiz),q21(mvsiz),
125 . f,df,yld_i,h(mvsiz),
126 . fisokin,hk(mvsiz),fhk,fa01,fa02,fa03,nu1,xfac,yfac,ce,
127 . e1(mvsiz), a11(mvsiz),a21(mvsiz),g1(mvsiz),g31(mvsiz),
128 . dydxe(mvsiz),escale(mvsiz),einf,normxx,normyy
129 my_real
130 . temp(mvsiz), vol0, rhocp
131 my_real,
132 . DIMENSION(NEL,3) :: xx
133 INTEGER,DIMENSION(NEL,3) :: IPOS
134C
135 DATA NMAX/2/
136C-----------------------------------------------
137C USER VARIABLES INITIALIZATION
138C-----------------------------------------------
139 FUNC_TAB = itable(1)
140 fisokin = uparam(1)
141 opte = uparam(23)
142 ce = uparam(25)
143c--------Two options with or without young evolution with
144c--------plastic strain
145 IF (opte==1.OR. ce > zero)THEN
146 DO i=1,nel
147 e1(i) = uparam(2)
148 a11(i) = uparam(3)
149 a21(i) = uparam(4)
150 g1(i) = uparam(5)
151 g31(i) = uparam(17)
152 ENDDO
153 nu = uparam(6)
154 a01 = uparam(7)
155 a02 = uparam(8)
156 a03 = uparam(9)
157 a12 = uparam(10)
158 epsmax = uparam(13)
159 einf = uparam(24)
160 IF(epsmax==zero)THEN
161c NXK=SIZE(TABLE(ITABLE)%X(1)%VALUES)
162c IF(TABLE(ITABLE)%Y%VALUES(NXK)==ZERO)THEN
163c EPSMAX = TABLE(ITABLE)%X(1)%VALUES(NXK)
164c ELSE
165 epsmax = ep30
166c ENDIF
167 ENDIF
168C
169 IF (opte == 1)THEN
170 DO i=1,nel
171 IF(pla(i) > zero)THEN
172 ifunce = uparam(22)
173 escale(i) = finter(ifunc(ifunce),pla(i),npf,tf,dydxe(i))
174 e1(i) = escale(i)* e1(i)
175 a11(i) = e1(i)/(one - nu*nu)
176 a21(i) = nu*a11(i)
177 g1(i) = half*e1(i)/(one+nu)
178 g31(i) = three*g1(i)
179 gs(i) = g1(i) * shf(i)
180 ENDIF
181 ENDDO
182 ELSEIF ( ce /= zero) THEN
183 DO i=1,nel
184 IF(pla(i) > zero)THEN
185 e1(i) = e1(i)-(e1(i)-einf)*(one-exp(-ce*pla(i)))
186 a11(i) = e1(i)/(one - nu*nu)
187 a21(i) = nu*a11(i)
188 g1(i) = half*e1(i)/(one+nu)
189 g31(i) = three*g1(i)
190 gs(i) = g1(i) * shf(i)
191 ENDIF
192 ENDDO
193 ENDIF
194 xfac =uparam(11)
195 yfac =uparam(12)
196C
197 nnu1 = nu / (one - nu)
198 epsr1 =uparam(14)
199 epsr2 =uparam(15)
200C
201 DO i=1,nel
202! SIGOXX(I) = SIGOXX(I) - UVAR(I,2)
203! SIGOYY(I) = SIGOYY(I) - UVAR(I,3)
204! SIGOXY(I) = SIGOXY(I) - UVAR(I,4)
205 dpla1(i) = zero
206 ENDDO
207C
208C-----------------------------------------------
209C calculation of the temperature (conduction or adiabatic)
210C--------------------
211 IF(jthe > 0 ) THEN
212 DO i=1,nel
213 temp(i) = tempel(i)
214 ENDDO
215 ELSE
216 rhocp=uparam(21)
217 DO i=1,nel
218 temp(i) = uparam(20)
219 vol0 = vol(i) * rho0(i)
220 temp(i) = temp(i) + (eint(i,1)+ eint(i,2)) * rhocp/vol0
221 ENDDO
222 ENDIF
223C
224 DO i=1,nel
225 signxx(i)=sigoxx(i) - uvar(i,2) +a11(i)*depsxx(i)+a21(i)*depsyy(i) !
226 signyy(i)=sigoyy(i) - uvar(i,3) +a21(i)*depsxx(i)+a11(i)*depsyy(i)!
227 signxy(i)=sigoxy(i) - uvar(i,4) +g1(i) *depsxy(i)!
228 signyz(i)=sigoyz(i)+gs(i)*depsyz(i)
229 signzx(i)=sigozx(i)+gs(i)*depszx(i)
230C
231 soundsp(i) = sqrt(a11(i)/rho0(i))!
232 viscmax(i) = zero
233 etse(i) = one
234C-------------------
235C STRAIN RATE
236C-------------------
237 epsp(i) = half*( abs(epspxx(i)+epspyy(i))
238 . + sqrt( (epspxx(i)-epspyy(i))*(epspxx(i)-epspyy(i))
239 . + epspxy(i)*epspxy(i) ) )
240C-------------------
241C STRAIN
242C-------------------
243 epst = half*( epsxx(i)+epsyy(i)
244 . + sqrt( (epsxx(i)-epsyy(i))*(epsxx(i)-epsyy(i))
245 . + epsxy(i)*epsxy(i) ) )
246 fail(i) = max(zero,min(one,(epsr2-epst)/(epsr2-epsr1)))
247C
248 ENDDO
249C-------------------
250C HARDENING LAW
251C-------------------
252 DO i=1,nel
253 ipos(i,1) = nint(uvar(i,5))
254 ipos(i,2) = nint(uvar(i,6))
255 ipos(i,3) = nint(uvar(i,7))
256C
257 xx(i,1)=pla(i)
258 xx(i,2)=epsp(i)*xfac
259 xx(i,3)=temp(i)
260 END DO
261C
262 CALL table_vinterp(table(func_tab),nel,nel,ipos,xx,yld,dydx)
263C
264
265 DO i=1,nel
266 yld(i) = yfac*yld(i)
267 yld(i) = fail(i)*yld(i)
268 yld(i) = max(yld(i),em20)
269 h(i) = fail(i)*dydx(i)
270 uvar(i,5) = ipos(i,1)
271 uvar(i,6) = ipos(i,2)
272 uvar(i,7) = ipos(i,3)
273 END DO
274C-----
275 IF(fisokin/=zero)THEN
276 DO i=1,nel
277 ipos(i,1) = 1
278C
279 xx(i,1)=zero
280 xx(i,2)=epsp(i)*xfac
281 xx(i,3)=temp(i)
282 END DO
283 CALL table_vinterp(table(func_tab),nel,nel,ipos,xx,yk,dydx)
284
285 DO i=1,nel
286C ECROUISSAGE CINEMATIQUE
287 yld(i) = (one-fisokin) * yld(i) +
288 . fisokin * fail(i) * yfac * yk(i)
289 yld(i) = max(yld(i),em20)
290 ENDDO
291 END IF
292C-------------------------
293C HILL CRITERION
294C-------------------------
295 DO i=1,nel
296 sigy(i) = yld(i)
297 h(i) = max(zero,h(i))
298 s1=a01*signxx(i)*signxx(i)
299 s2=a02*signyy(i)*signyy(i)
300 s3=a03*signxx(i)*signyy(i)
301 axy(i)=a12*signxy(i)*signxy(i)
302 svm(i)=sqrt(s1+s2-s3+axy(i))
303 IF (inloc == 0) THEN
304 dezz = -(depsxx(i)+depsyy(i))*nnu1
305 thk(i) = thk(i) + dezz*thkly(i)*off(i)
306 ENDIF
307 ENDDO
308C-------------------------
309C GATHER PLASTIC FLOW
310C-------------------------
311 nindx=0
312 DO i=1,nel
313 IF(svm(i)>yld(i).AND.off(i)==one) THEN
314 nindx=nindx+1
315 index(nindx)=i
316 ENDIF
317 ENDDO
318 IF(nindx>0) THEN
319C---------------------------
320C DEP EN CONTRAINTE PLANE
321C---------------------------
322 DO j=1,nindx
323 i=index(j)
324 dpla_j(i)=(svm(i)-yld(i))/(g31(i)+h(i)) !!!
325 etse(i)= h(i)/(h(i)+e1(i)) !!!
326C +++ This part can be Neglige when FHK is very small
327 hk(i) = h(i)*fisokin
328 fhk = four_over_3*hk(i)/a11(i) !!
329 fa01 =a01*fhk
330 fa02 =a02*fhk
331 fa03 =a03*fhk
332 nu1 =nu+half*fhk
333C ---
334 nu2 = 1.-nu1*nu1+fhk*fhk
335 nu3 = nu1*half
336 nu4 = half*(one-nu)
337 nu5 = one-nnu1
338 s1=a01*nu1*two-a03-fa03
339 s2=a02*nu1*two-a03-fa03
340 s12=a03-nu1*(a01+a02)+fa03
341 s3=sqrt(nu2*(a01-a02)*(a01-a02)+s12*s12)
342 IF (abs(s1)<em20) THEN
343 q12(i)=zero
344 ELSE
345 q12(i)=-(a01-a02+s3+fa01-fa02)/s1
346 ENDIF
347 IF (abs(s2)<em20) THEN
348 q21(i)=zero
349 ELSE
350 q21(i)=(a01-a02+s3+fa01-fa02)/s2
351 ENDIF
352 jq(i)=one/(1-q12(i)*q21(i))
353 jq2(i)=jq(i)*jq(i)
354 a=a01*q12(i)
355 b=a02*q21(i)
356 a_1=(a01+a03*q21(i)+b*q21(i))*jq2(i)
357 a_2=(a02+a03*q12(i)+a*q12(i))*jq2(i)
358 a_3=(a+b)*jq2(i)*two+a03*(jq2(i)*two-jq(i))
359 s11(i)=signxx(i)+signyy(i)*q12(i)
360 s22(i)=q21(i)*signxx(i)+signyy(i)
361 axx(i)=a_1*s11(i)*s11(i)
362 ayy(i)=a_2*s22(i)*s22(i)
363 a_xy(i)=a_3*s11(i)*s22(i)
364 a=a03*nu3
365 b=s3*jq(i)
366 b_1(i)=a02-a-b+fa02
367 b_2(i)=a01-a+b+fa01
368 b_3(i)=a12*(nu4+half*fhk)
369 h(i) = h(i)-hk(i)
370 h(i) = max(zero,h(i))
371 ENDDO
372C
373 DO n=1,nmax
374 DO j=1,nindx
375 i=index(j)
376 dpla_i(i)=dpla_j(i)
377 IF(dpla_i(i)>zero) THEN
378 yld_i =yld(i)+h(i)*dpla_i(i)
379 dr =a11(i)*dpla_i(i)/yld_i !
380 p_1=one/(one + b_1(i)*dr)
381 pp1=p_1*p_1
382 p_2=one/(one+b_2(i)*dr)
383 pp2=p_2*p_2
384 p_3=one/(one+b_3(i)*dr)
385 pp3=p_3*p_3
386 f =axx(i)*pp1+ayy(i)*pp2-a_xy(i)*p_1*p_2+axy(i)*pp3
387 . -yld_i*yld_i
388 df =-((axx(i)*p_1-a_xy(i)*p_2*half)*pp1*b_1(i)+
389 . (ayy(i)*p_2-a_xy(i)*p_1*half)*pp2*b_2(i)+
390 . axy(i)*pp3*p_3*b_3(i))*(a11(i)-dr*h(i))/yld_i
391 . -h(i)*yld_i
392 dpla_j(i)=max(zero,dpla_i(i)-f*half/df)
393 ELSE
394 dpla_j(i)=zero
395 ENDIF
396 ENDDO
397 ENDDO
398
399C------------------------------------------
400C CONTRAINTES PLASTIQUEMENT ADMISSIBLES
401C------------------------------------------
402 DO j=1,nindx
403 i=index(j)
404 pla(i) = pla(i) + dpla_j(i)
405 yld_i =yld(i)+h(i)*dpla_j(i)
406 dr0 =dpla_j(i)/yld_i
407 dr =a11(i)*dr0 !!!!
408 p_1=one/(one + b_1(i)*dr)
409 p_2=one/(one + b_2(i)*dr)
410 p_3=one/(one + b_3(i)*dr)
411 s1=s11(i)*p_1
412 s2=s22(i)*p_2
413 signxx(i)=jq(i)*(s1-s2*q12(i))
414 signyy(i)=jq(i)*(s2-s1*q21(i))
415 signxy(i)=signxy(i)*p_3
416 s1=a01*signxx(i)+a02*signyy(i)
417 . -a03*(signxx(i)+signyy(i))*half
418 IF (inloc == 0) THEN
419 dezz = - nu5*dpla_j(i)*s1/yld_i
420 thk(i) = thk(i) + dezz*thkly(i)*off(i)
421 ENDIF
422 s1 =a03*half
423 p_1=a01*signxx(i)-s1*signyy(i)
424 p_2=a02*signyy(i)-s1*signxx(i)
425 p_3=a12*signxy(i)
426 dr0 = two_third*dr0*hk(i)
427 uvar(i,2) = uvar(i,2) + (two*p_1+p_2)*dr0
428 uvar(i,3) = uvar(i,3) + (two*p_2+p_1)*dr0
429 uvar(i,4) = uvar(i,4) + half*p_3*dr0
430 dpla1(i) = dpla_j(i)
431 ENDDO
432C
433 DO i=1,nel
434 IF(pla(i)>epsmax.AND.off(i)==one)off(i)=four_over_5
435 ENDDO
436 ENDIF
437 DO i=1,nel
438 signxx(i) = signxx(i) + uvar(i,2)
439 signyy(i) = signyy(i) + uvar(i,3)
440 signxy(i) = signxy(i) + uvar(i,4)
441 ENDDO
442
443 ELSE
444 e = uparam(2)
445 a1 = uparam(3)
446 a2 = uparam(4)
447 g = uparam(5)
448 nu = uparam(6)
449 a01 = uparam(7)
450 a02 = uparam(8)
451 a03 = uparam(9)
452 a12 = uparam(10)
453 g3 = uparam(17)
454 epsmax = uparam(13)
455 IF(epsmax==zero)THEN
456c NXK=SIZE(TABLE(ITABLE)%X(1)%VALUES)
457c IF(TABLE(ITABLE)%Y%VALUES(NXK)==ZERO)THEN
458c EPSMAX = TABLE(ITABLE)%X(1)%VALUES(NXK)
459c ELSE
460 epsmax = ep30
461c ENDIF
462 ENDIF
463C
464 xfac =uparam(11)
465 yfac =uparam(12)
466C
467 nnu1 = nu / (one - nu)
468 epsr1 =uparam(14)
469 epsr2 =uparam(15)
470C
471 DO i=1,nel
472! SIGOXX(I) = SIGOXX(I) - UVAR(I,2)
473! SIGOYY(I) = SIGOYY(I) - UVAR(I,3)
474! SIGOXY(I) = SIGOXY(I) - UVAR(I,4)
475 dpla1(i) = zero
476 ENDDO
477C
478C-----------------------------------------------
479C calculation of the temperature (conduction or adiabatic)
480C--------------------
481 IF(jthe > 0 ) THEN
482 DO i=1,nel
483 temp(i) = tempel(i)
484 ENDDO
485 ELSE
486 rhocp=uparam(21)
487 DO i=1,nel
488 temp(i) = uparam(20)
489 vol0 = vol(i) * rho0(i)
490 temp(i) = temp(i) + (eint(i,1)+ eint(i,2)) * rhocp/vol0
491 ENDDO
492 ENDIF
493C
494 DO i=1,nel
495 signxx(i)=sigoxx(i) - uvar(i,2) +a1*depsxx(i)+a2*depsyy(i)
496 signyy(i)=sigoyy(i) - uvar(i,3) +a2*depsxx(i)+a1*depsyy(i)
497 signxy(i)=sigoxy(i) - uvar(i,4) +g *depsxy(i)
498 signyz(i)=sigoyz(i)+gs(i)*depsyz(i)
499 signzx(i)=sigozx(i)+gs(i)*depszx(i)
500C
501 soundsp(i) = sqrt(a1/rho0(i))
502 viscmax(i) = zero
503 etse(i) = one
504C-------------------
505C STRAIN RATE
506C-------------------
507 epsp(i) = half*( abs(epspxx(i)+epspyy(i))
508 . + sqrt( (epspxx(i)-epspyy(i))*(epspxx(i)-epspyy(i))
509 . + epspxy(i)*epspxy(i) ) )
510C-------------------
511C STRAIN
512C-------------------
513 epst = half*( epsxx(i)+epsyy(i)
514 . + sqrt( (epsxx(i)-epsyy(i))*(epsxx(i)-epsyy(i))
515 . + epsxy(i)*epsxy(i) ) )
516 fail(i) = max(zero,min(one,(epsr2-epst)/(epsr2-epsr1)))
517C
518 ENDDO
519C-------------------
520C HARDENING LAW
521C-------------------
522 DO i=1,nel
523 ipos(i,1) = nint(uvar(i,5))
524 ipos(i,2) = nint(uvar(i,6))
525 ipos(i,3) = nint(uvar(i,7))
526C
527 xx(i,1)=pla(i)
528 xx(i,2)=epsp(i)*xfac
529 xx(i,3)=temp(i)
530 END DO
531C
532 CALL table_vinterp(table(func_tab),nel,nel,ipos,xx,yld,dydx)
533C
534 DO i=1,nel
535 yld(i) = yfac*yld(i)
536 yld(i) = fail(i)*yld(i)
537 yld(i) = max(yld(i),em20)
538 h(i) = fail(i)*dydx(i)
539 uvar(i,5) = ipos(i,1)
540 uvar(i,6) = ipos(i,2)
541 uvar(i,7) = ipos(i,3)
542 END DO
543C-----
544 IF(fisokin/=zero)THEN
545 DO i=1,nel
546 ipos(i,1) = 1
547C
548 xx(i,1)=zero
549 xx(i,2)=epsp(i)*xfac
550 xx(i,3)=temp(i)
551 END DO
552 CALL table_vinterp(table(func_tab),nel,nel,ipos,xx,yk,dydx)
553
554 DO i=1,nel
555C ECROUISSAGE CINEMATIQUE
556 yld(i) = (one-fisokin) * yld(i) +
557 . fisokin * fail(i) * yfac * yk(i)
558 yld(i) = max(yld(i),em20)
559 ENDDO
560 END IF
561C-------------------------
562C HILL CRITERION
563C-------------------------
564 DO i=1,nel
565 sigy(i) = yld(i)
566 h(i) = max(zero,h(i))
567 s1=a01*signxx(i)*signxx(i)
568 s2=a02*signyy(i)*signyy(i)
569 s3=a03*signxx(i)*signyy(i)
570 axy(i)=a12*signxy(i)*signxy(i)
571 svm(i)=sqrt(s1+s2-s3+axy(i))
572 IF (inloc == 0) THEN
573 dezz = -(depsxx(i)+depsyy(i))*nnu1
574 thk(i) = thk(i) + dezz*thkly(i)*off(i)
575 ENDIF
576 ENDDO
577C-------------------------
578C GATHER PLASTIC FLOW
579C-------------------------
580 nindx=0
581 DO i=1,nel
582 IF(svm(i)>yld(i).AND.off(i)==one) THEN
583 nindx=nindx+1
584 index(nindx)=i
585 ENDIF
586 ENDDO
587 IF(nindx>0) THEN
588C---------------------------
589C DEP EN CONTRAINTE PLANE
590C---------------------------
591 DO j=1,nindx
592 i=index(j)
593 dpla_j(i)=(svm(i)-yld(i))/(g3+h(i))
594 etse(i)= h(i)/(h(i)+e)
595C +++ This part can be Neglige when FHK is very small
596 hk(i) = h(i)*fisokin
597 fhk = four_over_3*hk(i)/a1
598 fa01 =a01*fhk
599 fa02 =a02*fhk
600 fa03 =a03*fhk
601 nu1 =nu+half*fhk
602C ---
603 nu2 = 1.-nu1*nu1+fhk*fhk
604 nu3 = nu1*half
605 nu4 = half*(one-nu)
606 nu5 = one-nnu1
607 s1=a01*nu1*two-a03-fa03
608 s2=a02*nu1*two-a03-fa03
609 s12=a03-nu1*(a01+a02)+fa03
610 s3=sqrt(nu2*(a01-a02)*(a01-a02)+s12*s12)
611 IF (abs(s1)<em20) THEN
612 q12(i)=zero
613 ELSE
614 q12(i)=-(a01-a02+s3+fa01-fa02)/s1
615 ENDIF
616 IF (abs(s2)<em20) THEN
617 q21(i)=zero
618 ELSE
619 q21(i)=(a01-a02+s3+fa01-fa02)/s2
620 ENDIF
621 jq(i)=one/(1-q12(i)*q21(i))
622 jq2(i)=jq(i)*jq(i)
623 a=a01*q12(i)
624 b=a02*q21(i)
625 a_1=(a01+a03*q21(i)+b*q21(i))*jq2(i)
626 a_2=(a02+a03*q12(i)+a*q12(i))*jq2(i)
627 a_3=(a+b)*jq2(i)*two+a03*(jq2(i)*two-jq(i))
628 s11(i)=signxx(i)+signyy(i)*q12(i)
629 s22(i)=q21(i)*signxx(i)+signyy(i)
630 axx(i)=a_1*s11(i)*s11(i)
631 ayy(i)=a_2*s22(i)*s22(i)
632 a_xy(i)=a_3*s11(i)*s22(i)
633 a=a03*nu3
634 b=s3*jq(i)
635 b_1(i)=a02-a-b+fa02
636 b_2(i)=a01-a+b+fa01
637 b_3(i)=a12*(nu4+half*fhk)
638 h(i) = h(i)-hk(i)
639 h(i) = max(zero,h(i))
640 ENDDO
641C
642 DO n=1,nmax
643 DO j=1,nindx
644 i=index(j)
645 dpla_i(i)=dpla_j(i)
646 IF(dpla_i(i)>zero) THEN
647 yld_i =yld(i)+h(i)*dpla_i(i)
648 dr =a1*dpla_i(i)/yld_i
649 p_1=one/(one + b_1(i)*dr)
650 pp1=p_1*p_1
651 p_2=one/(one+b_2(i)*dr)
652 pp2=p_2*p_2
653 p_3=one/(one+b_3(i)*dr)
654 pp3=p_3*p_3
655 f =axx(i)*pp1+ayy(i)*pp2-a_xy(i)*p_1*p_2+axy(i)*pp3
656 . -yld_i*yld_i
657 df =-((axx(i)*p_1-a_xy(i)*p_2*half)*pp1*b_1(i)+
658 . (ayy(i)*p_2-a_xy(i)*p_1*half)*pp2*b_2(i)+
659 . axy(i)*pp3*p_3*b_3(i))*(a1-dr*h(i))/yld_i
660 . -h(i)*yld_i
661 dpla_j(i)=max(zero,dpla_i(i)-f*half/df)
662 ELSE
663 dpla_j(i)=zero
664 ENDIF
665 ENDDO
666 ENDDO
667C------------------------------------------
668C CONTRAINTES PLASTIQUEMENT ADMISSIBLES
669C------------------------------------------
670 DO j=1,nindx
671 i=index(j)
672 pla(i) = pla(i) + dpla_j(i)
673 yld_i =yld(i)+h(i)*dpla_j(i)
674 dr0 =dpla_j(i)/yld_i
675 dr =a1*dr0
676 p_1=one/(one + b_1(i)*dr)
677 p_2=one/(one + b_2(i)*dr)
678 p_3=one/(one + b_3(i)*dr)
679 s1=s11(i)*p_1
680 s2=s22(i)*p_2
681 signxx(i)=jq(i)*(s1-s2*q12(i))
682 signyy(i)=jq(i)*(s2-s1*q21(i))
683 signxy(i)=signxy(i)*p_3
684 s1=a01*signxx(i)+a02*signyy(i)
685 . -a03*(signxx(i)+signyy(i))*half
686 IF (inloc == 0) THEN
687 dezz = - nu5*dpla_j(i)*s1/yld_i
688 thk(i) = thk(i) + dezz*thkly(i)*off(i)
689 ENDIF
690 s1 =a03*half
691 p_1=a01*signxx(i)-s1*signyy(i)
692 p_2=a02*signyy(i)-s1*signxx(i)
693 p_3=a12*signxy(i)
694 dr0 = two_third*dr0*hk(i)
695 uvar(i,2) = uvar(i,2) + (two*p_1+p_2)*dr0
696 uvar(i,3) = uvar(i,3) + (two*p_2+p_1)*dr0
697 uvar(i,4) = uvar(i,4) + half*p_3*dr0
698 dpla1(i) = dpla_j(i)
699 ENDDO
700C
701 DO i=1,nel
702 IF(pla(i)>epsmax.AND.off(i)==one)off(i)=four_over_5
703 ENDDO
704 ENDIF
705 DO i=1,nel
706 signxx(i) = signxx(i) + uvar(i,2)
707 signyy(i) = signyy(i) + uvar(i,3)
708 signxy(i) = signxy(i) + uvar(i,4)
709 ENDDO
710 ENDIF ! (opte=1 or ce>0)
711C
712C--------------------------------
713C HARDENING MODULUS
714C--------------------------------
715 DO i=1,nel
716 hardm(i) = h(i)
717 ENDDO
718C------------------------------------------------------------------
719C EQUIVALENT STRESS FOR OUTPUT - NON-LOCAL THICKNESS VARIATION
720C------------------------------------------------------------------
721 DO i = 1,nel
722 s1 = a01*signxx(i)*signxx(i)
723 s2 = a02*signyy(i)*signyy(i)
724 s3 = a03*signxx(i)*signyy(i)
725 axy(i) = a12*signxy(i)*signxy(i)
726 svm(i) = sqrt(s1+s2-s3+axy(i))
727 seq_output(i) = svm(i)
728 IF ((inloc > 0).AND.(loff(i) == one)) THEN
729 normxx = (two*a01*signxx(i) - a03*signyy(i))/(max(two*svm(i),em20))
730 normyy = (two*a02*signyy(i) - a03*signxx(i))/(max(two*svm(i),em20))
731 dezz = max(dplanl(i),zero)*(normxx + normyy)
732 dezz = -nu*((signxx(i)-sigoxx(i)+signyy(i)-sigoyy(i))/e) - dezz
733 thk(i) = thk(i) + dezz*thkly(i)*off(i)
734 ENDIF
735 ENDDO
736!-------------------------
737 RETURN
738 END
area
subroutine area(d1, x, x2, y, y2, eint, stif0)
Definition hm_read_prop32.F:567
interface_table_mod::table_vinterp
Definition table_mod.F:252
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:239
table_mod
Definition table_mod.F:113
sigeps73c
subroutine sigeps73c(nel, nuparam, nuvar, 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, itable, etse, gs, sigy, dpla1, epsp, table, vol, tempel, die, npf, nfunc, ifunc, tf, shf, hardm, seq_output, inloc, dplanl, jthe, loff)
Definition sigeps73c.F:53
table_mod::ttable
Definition table_mod.F:126