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, COEF, 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),coef(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,NXK,J1,J2,N,NINDX,NMAX,
111 . INDEX(MVSIZ),OPTE,IFUNCE,FUNC_TAB
112 my_real
113 . A,B,S12,FAC,DEZZ,SIGZ,S1,S2,S3,DPLA,VM2,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 . svm0(mvsiz),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 . err,f,df,pla_i,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 Calcul de la tempurature. (conduction ou adiabatique)
210C--------------------
211 DO i=1,nel
212 coef(i) = one
213 END DO
214C
215 IF(jthe > 0 ) THEN
216 DO i=1,nel
217 temp(i) = tempel(i)
218 ENDDO
219 ELSE
220 rhocp=uparam(21)
221 DO i=1,nel
222 temp(i) = uparam(20)
223 vol0 = vol(i) * rho0(i)
224 temp(i) = temp(i)
225 . + coef(i)*rhocp*(eint(i,1)+ eint(i,2))/vol0
226 ENDDO
227 ENDIF
228C
229 DO i=1,nel
230 signxx(i)=sigoxx(i) - uvar(i,2) +a11(i)*depsxx(i)+a21(i)*depsyy(i) !
231 signyy(i)=sigoyy(i) - uvar(i,3) +a21(i)*depsxx(i)+a11(i)*depsyy(i)!
232 signxy(i)=sigoxy(i) - uvar(i,4) +g1(i) *depsxy(i)!
233 signyz(i)=sigoyz(i)+gs(i)*depsyz(i)
234 signzx(i)=sigozx(i)+gs(i)*depszx(i)
235C
236 soundsp(i) = sqrt(a11(i)/rho0(i))!
237 viscmax(i) = zero
238 etse(i) = one
239C-------------------
240C STRAIN RATE
241C-------------------
242 epsp(i) = half*( abs(epspxx(i)+epspyy(i))
243 . + sqrt( (epspxx(i)-epspyy(i))*(epspxx(i)-epspyy(i))
244 . + epspxy(i)*epspxy(i) ) )
245C-------------------
246C STRAIN
247C-------------------
248 epst = half*( epsxx(i)+epsyy(i)
249 . + sqrt( (epsxx(i)-epsyy(i))*(epsxx(i)-epsyy(i))
250 . + epsxy(i)*epsxy(i) ) )
251 fail(i) = max(zero,min(one,(epsr2-epst)/(epsr2-epsr1)))
252C
253 ENDDO
254C-------------------
255C HARDENING LAW
256C-------------------
257 DO i=1,nel
258 ipos(i,1) = nint(uvar(i,5))
259 ipos(i,2) = nint(uvar(i,6))
260 ipos(i,3) = nint(uvar(i,7))
261C
262 xx(i,1)=pla(i)
263 xx(i,2)=epsp(i)*xfac
264 xx(i,3)=temp(i)
265 END DO
266C
267 CALL table_vinterp(table(func_tab),nel,nel,ipos,xx,yld,dydx)
268C
269
270 DO i=1,nel
271 yld(i) = yfac*yld(i)
272 yld(i) = fail(i)*yld(i)
273 yld(i) = max(yld(i),em20)
274 h(i) = fail(i)*dydx(i)
275 uvar(i,5) = ipos(i,1)
276 uvar(i,6) = ipos(i,2)
277 uvar(i,7) = ipos(i,3)
278 END DO
279C-----
280 IF(fisokin/=zero)THEN
281 DO i=1,nel
282 ipos(i,1) = 1
283C
284 xx(i,1)=zero
285 xx(i,2)=epsp(i)*xfac
286 xx(i,3)=temp(i)
287 END DO
288 CALL table_vinterp(table(func_tab),nel,nel,ipos,xx,yk,dydx)
289
290 DO i=1,nel
291C ECROUISSAGE CINEMATIQUE
292 yld(i) = (one-fisokin) * yld(i) +
293 . fisokin * fail(i) * yfac * yk(i)
294 yld(i) = max(yld(i),em20)
295 ENDDO
296 END IF
297C-------------------------
298C HILL CRITERION
299C-------------------------
300 DO i=1,nel
301 sigy(i) = yld(i)
302 h(i) = max(zero,h(i))
303 s1=a01*signxx(i)*signxx(i)
304 s2=a02*signyy(i)*signyy(i)
305 s3=a03*signxx(i)*signyy(i)
306 axy(i)=a12*signxy(i)*signxy(i)
307 svm(i)=sqrt(s1+s2-s3+axy(i))
308 IF (inloc == 0) THEN
309 dezz = -(depsxx(i)+depsyy(i))*nnu1
310 thk(i) = thk(i) + dezz*thkly(i)*off(i)
311 ENDIF
312 ENDDO
313C-------------------------
314C GATHER PLASTIC FLOW
315C-------------------------
316 nindx=0
317 DO i=1,nel
318 IF(svm(i)>yld(i).AND.off(i)==one) THEN
319 nindx=nindx+1
320 index(nindx)=i
321 ENDIF
322 ENDDO
323 IF(nindx>0) THEN
324C---------------------------
325C DEP EN CONTRAINTE PLANE
326C---------------------------
327 DO j=1,nindx
328 i=index(j)
329 dpla_j(i)=(svm(i)-yld(i))/(g31(i)+h(i)) !!!
330 etse(i)= h(i)/(h(i)+e1(i)) !!!
331C +++cette partie peut etre neglige quand FHK est tres petit
332 hk(i) = h(i)*fisokin
333 fhk = four_over_3*hk(i)/a11(i) !!
334 fa01 =a01*fhk
335 fa02 =a02*fhk
336 fa03 =a03*fhk
337 nu1 =nu+half*fhk
338C ---
339 nu2 = 1.-nu1*nu1+fhk*fhk
340 nu3 = nu1*half
341 nu4 = half*(one-nu)
342 nu5 = one-nnu1
343 s1=a01*nu1*two-a03-fa03
344 s2=a02*nu1*two-a03-fa03
345 s12=a03-nu1*(a01+a02)+fa03
346 s3=sqrt(nu2*(a01-a02)*(a01-a02)+s12*s12)
347 IF (abs(s1)<em20) THEN
348 q12(i)=zero
349 ELSE
350 q12(i)=-(a01-a02+s3+fa01-fa02)/s1
351 ENDIF
352 IF (abs(s2)<em20) THEN
353 q21(i)=zero
354 ELSE
355 q21(i)=(a01-a02+s3+fa01-fa02)/s2
356 ENDIF
357 jq(i)=one/(1-q12(i)*q21(i))
358 jq2(i)=jq(i)*jq(i)
359 a=a01*q12(i)
360 b=a02*q21(i)
361 a_1=(a01+a03*q21(i)+b*q21(i))*jq2(i)
362 a_2=(a02+a03*q12(i)+a*q12(i))*jq2(i)
363 a_3=(a+b)*jq2(i)*two+a03*(jq2(i)*two-jq(i))
364 s11(i)=signxx(i)+signyy(i)*q12(i)
365 s22(i)=q21(i)*signxx(i)+signyy(i)
366 axx(i)=a_1*s11(i)*s11(i)
367 ayy(i)=a_2*s22(i)*s22(i)
368 a_xy(i)=a_3*s11(i)*s22(i)
369 a=a03*nu3
370 b=s3*jq(i)
371 b_1(i)=a02-a-b+fa02
372 b_2(i)=a01-a+b+fa01
373 b_3(i)=a12*(nu4+half*fhk)
374 h(i) = h(i)-hk(i)
375 h(i) = max(zero,h(i))
376 ENDDO
377C
378 DO n=1,nmax
379 DO j=1,nindx
380 i=index(j)
381 dpla_i(i)=dpla_j(i)
382 IF(dpla_i(i)>zero) THEN
383 yld_i =yld(i)+h(i)*dpla_i(i)
384 dr =a11(i)*dpla_i(i)/yld_i !
385 p_1=one/(one + b_1(i)*dr)
386 pp1=p_1*p_1
387 p_2=one/(one+b_2(i)*dr)
388 pp2=p_2*p_2
389 p_3=one/(one+b_3(i)*dr)
390 pp3=p_3*p_3
391 f =axx(i)*pp1+ayy(i)*pp2-a_xy(i)*p_1*p_2+axy(i)*pp3
392 . -yld_i*yld_i
393 df =-((axx(i)*p_1-a_xy(i)*p_2*half)*pp1*b_1(i)+
394 . (ayy(i)*p_2-a_xy(i)*p_1*half)*pp2*b_2(i)+
395 . axy(i)*pp3*p_3*b_3(i))*(a11(i)-dr*h(i))/yld_i
396 . -h(i)*yld_i
397 dpla_j(i)=max(zero,dpla_i(i)-f*half/df)
398 ELSE
399 dpla_j(i)=zero
400 ENDIF
401 ENDDO
402 ENDDO
403
404C------------------------------------------
405C CONTRAINTES PLASTIQUEMENT ADMISSIBLES
406C------------------------------------------
407 DO j=1,nindx
408 i=index(j)
409 pla(i) = pla(i) + dpla_j(i)
410 yld_i =yld(i)+h(i)*dpla_j(i)
411 dr0 =dpla_j(i)/yld_i
412 dr =a11(i)*dr0 !!!!
413 p_1=one/(one + b_1(i)*dr)
414 p_2=one/(one + b_2(i)*dr)
415 p_3=one/(one + b_3(i)*dr)
416 s1=s11(i)*p_1
417 s2=s22(i)*p_2
418 signxx(i)=jq(i)*(s1-s2*q12(i))
419 signyy(i)=jq(i)*(s2-s1*q21(i))
420 signxy(i)=signxy(i)*p_3
421 s1=a01*signxx(i)+a02*signyy(i)
422 . -a03*(signxx(i)+signyy(i))*half
423 IF (inloc == 0) THEN
424 dezz = - nu5*dpla_j(i)*s1/yld_i
425 thk(i) = thk(i) + dezz*thkly(i)*off(i)
426 ENDIF
427 s1 =a03*half
428 p_1=a01*signxx(i)-s1*signyy(i)
429 p_2=a02*signyy(i)-s1*signxx(i)
430 p_3=a12*signxy(i)
431 dr0 = two_third*dr0*hk(i)
432 uvar(i,2) = uvar(i,2) + (two*p_1+p_2)*dr0
433 uvar(i,3) = uvar(i,3) + (two*p_2+p_1)*dr0
434 uvar(i,4) = uvar(i,4) + half*p_3*dr0
435 dpla1(i) = dpla_j(i)
436 ENDDO
437C
438 DO i=1,nel
439 IF(pla(i)>epsmax.AND.off(i)==one)off(i)=four_over_5
440 ENDDO
441 ENDIF
442 DO i=1,nel
443 signxx(i) = signxx(i) + uvar(i,2)
444 signyy(i) = signyy(i) + uvar(i,3)
445 signxy(i) = signxy(i) + uvar(i,4)
446 ENDDO
447
448 ELSE
449 e = uparam(2)
450 a1 = uparam(3)
451 a2 = uparam(4)
452 g = uparam(5)
453 nu = uparam(6)
454 a01 = uparam(7)
455 a02 = uparam(8)
456 a03 = uparam(9)
457 a12 = uparam(10)
458 g3 = uparam(17)
459 epsmax = uparam(13)
460 IF(epsmax==zero)THEN
461c NXK=SIZE(TABLE(ITABLE)%X(1)%VALUES)
462c IF(TABLE(ITABLE)%Y%VALUES(NXK)==ZERO)THEN
463c EPSMAX = TABLE(ITABLE)%X(1)%VALUES(NXK)
464c ELSE
465 epsmax = ep30
466c ENDIF
467 ENDIF
468C
469 xfac =uparam(11)
470 yfac =uparam(12)
471C
472 nnu1 = nu / (one - nu)
473 epsr1 =uparam(14)
474 epsr2 =uparam(15)
475C
476 DO i=1,nel
477! SIGOXX(I) = SIGOXX(I) - UVAR(I,2)
478! SIGOYY(I) = SIGOYY(I) - UVAR(I,3)
479! SIGOXY(I) = SIGOXY(I) - UVAR(I,4)
480 dpla1(i) = zero
481 ENDDO
482C
483C-----------------------------------------------
484C Calcul de la tempurature. (conduction ou adiabatique)
485C--------------------
486 DO i=1,nel
487 coef(i) = one
488 END DO
489C
490 IF(jthe > 0 ) THEN
491 DO i=1,nel
492 temp(i) = tempel(i)
493 ENDDO
494 ELSE
495 rhocp=uparam(21)
496 DO i=1,nel
497 temp(i) = uparam(20)
498 vol0 = vol(i) * rho0(i)
499 temp(i) = temp(i)
500 . + coef(i)*rhocp*(eint(i,1)+ eint(i,2))/vol0
501 ENDDO
502 ENDIF
503C
504 DO i=1,nel
505 signxx(i)=sigoxx(i) - uvar(i,2) +a1*depsxx(i)+a2*depsyy(i)
506 signyy(i)=sigoyy(i) - uvar(i,3) +a2*depsxx(i)+a1*depsyy(i)
507 signxy(i)=sigoxy(i) - uvar(i,4) +g *depsxy(i)
508 signyz(i)=sigoyz(i)+gs(i)*depsyz(i)
509 signzx(i)=sigozx(i)+gs(i)*depszx(i)
510C
511 soundsp(i) = sqrt(a1/rho0(i))
512 viscmax(i) = zero
513 etse(i) = one
514C-------------------
515C STRAIN RATE
516C-------------------
517 epsp(i) = half*( abs(epspxx(i)+epspyy(i))
518 . + sqrt( (epspxx(i)-epspyy(i))*(epspxx(i)-epspyy(i))
519 . + epspxy(i)*epspxy(i) ) )
520C-------------------
521C STRAIN
522C-------------------
523 epst = half*( epsxx(i)+epsyy(i)
524 . + sqrt( (epsxx(i)-epsyy(i))*(epsxx(i)-epsyy(i))
525 . + epsxy(i)*epsxy(i) ) )
526 fail(i) = max(zero,min(one,(epsr2-epst)/(epsr2-epsr1)))
527C
528 ENDDO
529C-------------------
530C HARDENING LAW
531C-------------------
532 DO i=1,nel
533 ipos(i,1) = nint(uvar(i,5))
534 ipos(i,2) = nint(uvar(i,6))
535 ipos(i,3) = nint(uvar(i,7))
536C
537 xx(i,1)=pla(i)
538 xx(i,2)=epsp(i)*xfac
539c print *,xx(i,2)
540 xx(i,3)=temp(i)
541 END DO
542C
543 CALL table_vinterp(table(func_tab),nel,nel,ipos,xx,yld,dydx)
544C
545 DO i=1,nel
546 yld(i) = yfac*yld(i)
547 yld(i) = fail(i)*yld(i)
548 yld(i) = max(yld(i),em20)
549 h(i) = fail(i)*dydx(i)
550 uvar(i,5) = ipos(i,1)
551 uvar(i,6) = ipos(i,2)
552 uvar(i,7) = ipos(i,3)
553 END DO
554C-----
555 IF(fisokin/=zero)THEN
556 DO i=1,nel
557 ipos(i,1) = 1
558C
559 xx(i,1)=zero
560 xx(i,2)=epsp(i)*xfac
561 xx(i,3)=temp(i)
562 END DO
563 CALL table_vinterp(table(func_tab),nel,nel,ipos,xx,yk,dydx)
564
565 DO i=1,nel
566C ECROUISSAGE CINEMATIQUE
567 yld(i) = (one-fisokin) * yld(i) +
568 . fisokin * fail(i) * yfac * yk(i)
569 yld(i) = max(yld(i),em20)
570 ENDDO
571 END IF
572C-------------------------
573C HILL CRITERION
574C-------------------------
575 DO i=1,nel
576 sigy(i) = yld(i)
577 h(i) = max(zero,h(i))
578 s1=a01*signxx(i)*signxx(i)
579 s2=a02*signyy(i)*signyy(i)
580 s3=a03*signxx(i)*signyy(i)
581 axy(i)=a12*signxy(i)*signxy(i)
582 svm(i)=sqrt(s1+s2-s3+axy(i))
583 IF (inloc == 0) THEN
584 dezz = -(depsxx(i)+depsyy(i))*nnu1
585 thk(i) = thk(i) + dezz*thkly(i)*off(i)
586 ENDIF
587 ENDDO
588C-------------------------
589C GATHER PLASTIC FLOW
590C-------------------------
591 nindx=0
592 DO i=1,nel
593 IF(svm(i)>yld(i).AND.off(i)==one) THEN
594 nindx=nindx+1
595 index(nindx)=i
596 ENDIF
597 ENDDO
598 IF(nindx>0) THEN
599C---------------------------
600C DEP EN CONTRAINTE PLANE
601C---------------------------
602 DO j=1,nindx
603 i=index(j)
604 dpla_j(i)=(svm(i)-yld(i))/(g3+h(i))
605 etse(i)= h(i)/(h(i)+e)
606C +++cette partie peut etre neglige quand FHK est tres petit
607 hk(i) = h(i)*fisokin
608 fhk = four_over_3*hk(i)/a1
609 fa01 =a01*fhk
610 fa02 =a02*fhk
611 fa03 =a03*fhk
612 nu1 =nu+half*fhk
613C ---
614 nu2 = 1.-nu1*nu1+fhk*fhk
615 nu3 = nu1*half
616 nu4 = half*(one-nu)
617 nu5 = one-nnu1
618 s1=a01*nu1*two-a03-fa03
619 s2=a02*nu1*two-a03-fa03
620 s12=a03-nu1*(a01+a02)+fa03
621 s3=sqrt(nu2*(a01-a02)*(a01-a02)+s12*s12)
622 IF (abs(s1)<em20) THEN
623 q12(i)=zero
624 ELSE
625 q12(i)=-(a01-a02+s3+fa01-fa02)/s1
626 ENDIF
627 IF (abs(s2)<em20) THEN
628 q21(i)=zero
629 ELSE
630 q21(i)=(a01-a02+s3+fa01-fa02)/s2
631 ENDIF
632 jq(i)=one/(1-q12(i)*q21(i))
633 jq2(i)=jq(i)*jq(i)
634 a=a01*q12(i)
635 b=a02*q21(i)
636 a_1=(a01+a03*q21(i)+b*q21(i))*jq2(i)
637 a_2=(a02+a03*q12(i)+a*q12(i))*jq2(i)
638 a_3=(a+b)*jq2(i)*two+a03*(jq2(i)*two-jq(i))
639 s11(i)=signxx(i)+signyy(i)*q12(i)
640 s22(i)=q21(i)*signxx(i)+signyy(i)
641 axx(i)=a_1*s11(i)*s11(i)
642 ayy(i)=a_2*s22(i)*s22(i)
643 a_xy(i)=a_3*s11(i)*s22(i)
644 a=a03*nu3
645 b=s3*jq(i)
646 b_1(i)=a02-a-b+fa02
647 b_2(i)=a01-a+b+fa01
648 b_3(i)=a12*(nu4+half*fhk)
649 h(i) = h(i)-hk(i)
650 h(i) = max(zero,h(i))
651 ENDDO
652C
653 DO n=1,nmax
654 DO j=1,nindx
655 i=index(j)
656 dpla_i(i)=dpla_j(i)
657 IF(dpla_i(i)>zero) THEN
658 yld_i =yld(i)+h(i)*dpla_i(i)
659 dr =a1*dpla_i(i)/yld_i
660 p_1=one/(one + b_1(i)*dr)
661 pp1=p_1*p_1
662 p_2=one/(one+b_2(i)*dr)
663 pp2=p_2*p_2
664 p_3=one/(one+b_3(i)*dr)
665 pp3=p_3*p_3
666 f =axx(i)*pp1+ayy(i)*pp2-a_xy(i)*p_1*p_2+axy(i)*pp3
667 . -yld_i*yld_i
668 df =-((axx(i)*p_1-a_xy(i)*p_2*half)*pp1*b_1(i)+
669 . (ayy(i)*p_2-a_xy(i)*p_1*half)*pp2*b_2(i)+
670 . axy(i)*pp3*p_3*b_3(i))*(a1-dr*h(i))/yld_i
671 . -h(i)*yld_i
672 dpla_j(i)=max(zero,dpla_i(i)-f*half/df)
673 ELSE
674 dpla_j(i)=zero
675 ENDIF
676 ENDDO
677 ENDDO
678C------------------------------------------
679C CONTRAINTES PLASTIQUEMENT ADMISSIBLES
680C------------------------------------------
681 DO j=1,nindx
682 i=index(j)
683 pla(i) = pla(i) + dpla_j(i)
684 yld_i =yld(i)+h(i)*dpla_j(i)
685 dr0 =dpla_j(i)/yld_i
686 dr =a1*dr0
687 p_1=one/(one + b_1(i)*dr)
688 p_2=one/(one + b_2(i)*dr)
689 p_3=one/(one + b_3(i)*dr)
690 s1=s11(i)*p_1
691 s2=s22(i)*p_2
692 signxx(i)=jq(i)*(s1-s2*q12(i))
693 signyy(i)=jq(i)*(s2-s1*q21(i))
694 signxy(i)=signxy(i)*p_3
695 s1=a01*signxx(i)+a02*signyy(i)
696 . -a03*(signxx(i)+signyy(i))*half
697 IF (inloc == 0) THEN
698 dezz = - nu5*dpla_j(i)*s1/yld_i
699 thk(i) = thk(i) + dezz*thkly(i)*off(i)
700 ENDIF
701 s1 =a03*half
702 p_1=a01*signxx(i)-s1*signyy(i)
703 p_2=a02*signyy(i)-s1*signxx(i)
704 p_3=a12*signxy(i)
705 dr0 = two_third*dr0*hk(i)
706 uvar(i,2) = uvar(i,2) + (two*p_1+p_2)*dr0
707 uvar(i,3) = uvar(i,3) + (two*p_2+p_1)*dr0
708 uvar(i,4) = uvar(i,4) + half*p_3*dr0
709 dpla1(i) = dpla_j(i)
710 ENDDO
711C
712 DO i=1,nel
713 IF(pla(i)>epsmax.AND.off(i)==one)off(i)=four_over_5
714 ENDDO
715 ENDIF
716 DO i=1,nel
717 signxx(i) = signxx(i) + uvar(i,2)
718 signyy(i) = signyy(i) + uvar(i,3)
719 signxy(i) = signxy(i) + uvar(i,4)
720 ENDDO
721 ENDIF ! (opte=1 or ce>0)
722C
723C--------------------------------
724C HARDENING MODULUS
725C--------------------------------
726 DO i=1,nel
727 hardm(i) = h(i)
728 ENDDO
729C------------------------------------------------------------------
730C EQUIVALENT STRESS FOR OUTPUT - NON-LOCAL THICKNESS VARIATION
731C------------------------------------------------------------------
732 DO i = 1,nel
733 s1 = a01*signxx(i)*signxx(i)
734 s2 = a02*signyy(i)*signyy(i)
735 s3 = a03*signxx(i)*signyy(i)
736 axy(i) = a12*signxy(i)*signxy(i)
737 svm(i) = sqrt(s1+s2-s3+axy(i))
738 seq_output(i) = svm(i)
739 IF ((inloc > 0).AND.(loff(i) == one)) THEN
740 normxx = (two*a01*signxx(i) - a03*signyy(i))/(max(two*svm(i),em20))
741 normyy = (two*a02*signyy(i) - a03*signxx(i))/(max(two*svm(i),em20))
742 dezz = max(dplanl(i),zero)*(normxx + normyy)
743 dezz = -nu*((signxx(i)-sigoxx(i)+signyy(i)-sigoyy(i))/e) - dezz
744 thk(i) = thk(i) + dezz*thkly(i)*off(i)
745 ENDIF
746 ENDDO
747!-------------------------
748 RETURN
749 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:251
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
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, coef, npf, nfunc, ifunc, tf, shf, hardm, seq_output, inloc, dplanl, jthe, loff)
Definition sigeps73c.F:53
table_mod::ttable
Definition table_mod.F:125