OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps38.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!|| sigeps38 ../engine/source/materials/mat/mat038/sigeps38.F
25!||--- called by ------------------------------------------------------
26!|| mulaw ../engine/source/materials/mat_share/mulaw.F90
27!|| mulaw8 ../engine/source/materials/mat_share/mulaw8.F90
28!||--- calls -----------------------------------------------------
29!|| checkaxes ../engine/source/materials/mat/mat038/sigeps38.F
30!|| dreh ../engine/source/materials/mat/mat033/sigeps33.F
31!|| finter ../engine/source/tools/curve/finter.F
32!|| jacobiew ../engine/source/materials/mat/mat033/sigeps33.F
33!|| valpvec_v ../engine/source/materials/mat/mat033/sigeps33.F
34!|| valpvecdp_v ../engine/source/materials/mat/mat033/sigeps33.F
35!||====================================================================
36 SUBROUTINE sigeps38(
37 1 NEL , NUPARAM, NUVAR , NFUNC , IFUNC , NPF ,
38 2 TF , TIME , TIMESTEP, UPARAM, RHO0 , RHO ,
39 3 VOLUME , EINT ,
40 4 EPSPXX , EPSPYY , EPSPZZ , EPSPXY, EPSPYZ, EPSPZX,
41 5 DEPSXX , DEPSYY , DEPSZZ , DEPSXY, DEPSYZ, DEPSZX,
42 6 EPSXX , EPSYY , EPSZZ , EPSXY , EPSYZ , EPSZX ,
43 7 SIGOXX , SIGOYY , SIGOZZ , SIGOXY, SIGOYZ, SIGOZX,
44 8 SIGNXX , SIGNYY , SIGNZZ , SIGNXY, SIGNYZ, SIGNZX,
45 9 SIGVXX , SIGVYY , SIGVZZ , SIGVXY, SIGVYZ, SIGVZX,
46 A SOUNDSP, VISCMAX, UVAR , OFF ,
47 B ISMSTR , MFXX , MFXY , MFXZ , MFYX ,
48 C MFYY , MFYZ , MFZX , MFZY , MFZZ , ET ,
49 D IHET ,OFFG ,EPSD )
50C-----------------------------------------------
51C I M P L I C I T T Y P E S
52C-----------------------------------------------
53#include "implicit_f.inc"
54C-----------------------------------------------
55C G l o b a l P a r a m e t e r s
56C-----------------------------------------------
57#include "mvsiz_p.inc"
58#include "scr05_c.inc"
59#include "impl1_c.inc"
60C----------------------------------------------------------------
61C I N P U T A R G U M E N T S
62C----------------------------------------------------------------
63 INTEGER NEL, NUPARAM, NUVAR,ISMSTR,IHET
64
65 my_real
66 . TIME , TIMESTEP , UPARAM(NUPARAM),
67 . RHO (NEL), RHO0 (NEL), VOLUME(NEL), EINT(NEL),
68 . EPSPXX(NEL), EPSPYY(NEL), EPSPZZ(NEL),
69 . EPSPXY(NEL), EPSPYZ(NEL), EPSPZX(NEL),
70 . DEPSXX(NEL), DEPSYY(NEL), DEPSZZ(NEL),
71 . DEPSXY(NEL), DEPSYZ(NEL), DEPSZX(NEL),
72 . EPSXX (NEL), EPSYY (NEL), EPSZZ (NEL),
73 . EPSXY (NEL), EPSYZ (NEL), EPSZX (NEL),
74 . SIGOXX(NEL), SIGOYY(NEL), SIGOZZ(NEL),
75 . SIGOXY(NEL), SIGOYZ(NEL), SIGOZX(NEL),
76 . mfxx(nel) , mfxy(nel), mfxz(nel),
77 . mfyx(nel) , mfyy(nel), mfyz(nel),
78 . mfzx(nel) , mfzy(nel), mfzz(nel),offg(nel)
79C----------------------------------------------------------------
80C O U T P U T A R G U M E N T S
81C----------------------------------------------------------------
82 my_real
83 . signxx(nel), signyy(nel), signzz(nel),
84 . signxy(nel), signyz(nel), signzx(nel),
85 . sigvxx(nel), sigvyy(nel), sigvzz(nel),
86 . sigvxy(nel), sigvyz(nel), sigvzx(nel),
87 . soundsp(nel), viscmax(nel),
88 . den1,den2,den3,den1d,den2d,
89 . den3d,old1,old2,old3,old4,old5,old6,et(*),epsd(nel)
90C----------------------------------------------------------------
91C I N P U T O U T P U T A R G U M E N T S
92C----------------------------------------------------------------
93 my_real
94 . uvar(nel,nuvar), off(nel)
95C----------------------------------------------------------------
96C VARIABLES FOR FUNCTION INTERPOLATION
97C----------------------------------------------------------------
98 INTEGER NPF(*), NFUNC, IFUNC(NFUNC)
99 my_real
100 . FINTER,TF(*)
101 EXTERNAL FINTER
102C----------------------------------------------------------------
103C L O C A L V A R I B L E S
104C----------------------------------------------------------------
105 INTEGER MFUNC,IUNLOAD
106 INTEGER NUPARAM0
107 INTEGER I,J,L,N,M,MDIR,NROT
108 INTEGER K(MVSIZ,3),K1(MVSIZ,3),KF(MVSIZ,3)
109 INTEGER KF1(MVSIZ,3),KUN(MVSIZ,3)
110 INTEGER IFLAG,ITOTAL,IMSTA
111 INTEGER NFUNC1,NFUNCUL,NFUNCP,NPCURVE,KCOMPAIR
112 INTEGER KRECOVER,KDECAY
113 LOGICAL TOTAL,UNCHANGED,JACOBI,UNLOADING
114 INTEGER IND1,LD
115 INTEGER II,JJ,KEN,IDEAC
116C REAL
117 my_real
118 . strain(mvsiz,3),rate(mvsiz,3),ratem(mvsiz),
119 . visc(mvsiz,3),
120c . ESECANT(MVSIZ,3),ETANGENT(MVSIZ,3),EELASTIC(MVSIZ,3),
121 . dsigma(mvsiz,3),decay,tensioncut,
122 . amin,amax,tolerance,lamda,efinal,epsfin,efactor,
123 . a(3,3),sn(3,3),earj(3),
124 . psc(mvsiz,3),ear(mvsiz,3),
125 . ebr(mvsiz,3),ecr(mvsiz,3),
126 . psn(mvsiz,3),ean(mvsiz,3),
127 . ebn(mvsiz,3),ecn(mvsiz,3),
128 . el(mvsiz,3),
129 . dpra(3,3),dprao(3,3),dprad(3,3),
130 . av(mvsiz,6),earv(mvsiz,3),dirprv(mvsiz,3,3),
131 . sigprv(mvsiz,3),ei(mvsiz),
132 . e0,vt,vc,rv,kkk,ggg,
133 . beta,hyster,ratedamp,theta,
134 . p0,relaxp,maxpres,phi,gamma,
135 . pair(mvsiz),volumer(mvsiz),
136 . df,df1,df2,
137 . viscosity,
138 . alpha(5),edot(5),
139 . psn1(mvsiz,3),psn2(mvsiz,3),
140 . edot0(mvsiz,3),edots(mvsiz,3),
141 . edotl(mvsiz,3),edotu(mvsiz,3),
142 . exponas,exponbs,funload,runload,
143 . e12(mvsiz),e23(mvsiz),e31(mvsiz),
144 . a11(mvsiz),a12(mvsiz),a13(mvsiz),a22(mvsiz),
145 . a23(mvsiz),a33(mvsiz),detc(mvsiz),
146 . v12(mvsiz),v23(mvsiz),v31(mvsiz),
147 . d11(mvsiz),d12(mvsiz),d44(mvsiz),
148 . emax(mvsiz),eyn(mvsiz,3),
149 . huge,small,tiny,shift,pscale
150 my_real
151 . ar1(6),dar1(6),ar2(6),dar2(6),dar3(3),earjd(3),
152 . tmp0,tmp1,tmp2,tmp3,tmp4,tmp5,tmp6,ad(3,3),
153 . aa(mvsiz),sigmax(mvsiz),esecant(mvsiz,3),
154 . tmp0p,tmp1p,tmp2p,tmp3p,tmp4p,tmp5p,tmp6p, dti,iemax,efac(mvsiz),
155 . ratio, pui,pui_1,pui_2,eymin,eymax,eyav
156C=======================================================================
157 huge = ep30
158 small= em3
159 tiny = em30
160C......................................................................
161C SET INITIAL MATERIAL CONSTANTS
162
163 nuparam0= uparam(1)
164
165 e0 = uparam(2)
166 vt = uparam(3)
167 vc = uparam(4)
168 rv = uparam(5)
169 iflag = uparam(6)
170 itotal = uparam(7)
171
172 beta = uparam(8)
173 hyster = uparam(9)
174 ratedamp = uparam(10)
175 krecover = uparam(11)
176 kdecay = uparam(12)
177 theta = uparam(13)
178c THETA = RATEDAMP
179
180 kcompair = uparam(14)
181 p0 = uparam(15)
182 gamma = uparam(16)
183 relaxp = uparam(17)
184 maxpres = uparam(18)
185 phi = uparam(19)
186
187 iunload = uparam(20)
188 funload = uparam(21)
189 runload = uparam(22)
190 exponas = uparam(23)
191 exponbs = uparam(24)
192
193 mfunc = uparam(25)
194 imsta = uparam(26)
195 ideac=0
196 IF (imsta>=10) THEN
197 imsta=imsta-10
198 ideac=1
199 END IF
200 tensioncut= uparam(27)
201
202 efinal = uparam(28)
203 epsfin = uparam(29)
204 lamda = uparam(30)
205 viscosity = uparam(31)
206 tolerance = uparam(32)
207 pscale = uparam(33)
208 nfunc1=(nfunc-2)/2
209* unloading function number
210 nfuncul=nfunc-1
211* function number for enclosed air pressure
212 nfuncp=nfunc
213
214* set strain rates and scale factors for curves
215 DO i=1,nfunc1
216 edot(i) =uparam(nuparam-nfunc1+i)
217 alpha(i)=uparam(nuparam-nfunc1*2+i)
218 ENDDO
219
220* initialize e-module
221 total=.true.
222C ITOTAL=0 ; TOTAL - LINEAR IN TENSION
223C ITOTAL=1 ; TOTAL - NONLINEAR IN TENSION
224C ITOTAL=2 ; INCREMENTAL - LINEAR IN TENSION
225C ITOTAL=3 ; INCREMENTAL - NONLINEAR IN TENSION
226 IF(itotal==2.OR.itotal==3) total=.false.
227 IF (ismstr>=10.AND.ismstr>=12) total=.true.
228CC +++ ([F]=[M_F]+[1]) ----- Done before
229C IF(ISMSTR==10) THEN
230CC--------- [B]=[F][F]^t strain-----
231C DO I=1,NEL
232C EPSXX(I)=MFXX(I)*(TWO+MFXX(I))+
233C . MFXY(I)*MFXY(I)+MFXZ(I)*MFXZ(I)
234C EPSYY(I)=MFYY(I)*(TWO+MFYY(I))+
235C . MFYX(I)*MFYX(I)+MFYZ(I)*MFYZ(I)
236C EPSZZ(I)=MFZZ(I)*(TWO+MFZZ(I))+
237C . MFZX(I)*MFZX(I)+MFZY(I)*MFZY(I)
238C EPSXY(I)=TWO*(MFXY(I)+MFYX(I)+MFXX(I)*MFYX(I)+
239C . MFXY(I)*MFYY(I)+MFXZ(I)*MFYZ(I))
240C EPSZX(I)=TWO*(MFXZ(I)+MFZX(I)+MFXX(I)*MFZX(I)+
241C . MFXY(I)*MFZY(I)+MFXZ(I)*MFZZ(I))
242C EPSYZ(I)=TWO*(MFZY(I)+MFYZ(I)+MFZX(I)*MFYX(I)+
243C . MFZY(I)*MFYY(I)+MFZZ(I)*MFYZ(I))
244C ENDDO
245C IF(IDTMIN(1)==3)THEN
246C DO I=1,NEL
247C IF (OFFG(I) <=ONE) CYCLE
248C EPSXX(I)=MFXX(I)
249C EPSYY(I)=MFYY(I)
250C EPSZZ(I)=MFZZ(I)
251C EPSXY(I)=MFXY(I)+MFYX(I)
252C EPSZX(I)=MFXZ(I)+MFZX(I)
253C EPSYZ(I)=MFZY(I)+MFYZ(I)
254C ENDDO
255C END IF
256C TOTAL=.TRUE.
257C ENDIF
258
259C......................................................................
260C INITIALIZE
261 DO i=1,nel
262 ei(i)=e0
263 DO j=1,3
264 visc(i,j)=zero
265 ENDDO
266C IF(TIME==ZERO) THEN
267C EYN(I,1) = E0
268C EYN(I,2) = E0
269C EYN(I,3) = E0
270C STRAIN(1-3),STRESS(4-6),STRAIN RATE(7-9)
271C DO J=1,9
272C UVAR(I,J)=ZERO
273C ENDDO
274C MODULE (EYN INITIALIZED TO E0)
275C DO J=10,12
276C UVAR(I,J)=E0
277C ENDDO
278C POISSON'S RATIO/MODULE
279C DO J=13,15
280C UVAR(I,J)=VT/E0
281C ENDDO
282C PRESSURE
283C UVAR(I,16)=ZERO
284C PRINCIPAL DIRECTIONS
285C DO J=17,25
286C UVAR(I,J)=ZERO
287C ENDDO
288C UVAR(I,17)=ONE
289C UVAR(I,21)=ONE
290C UVAR(I,25)=ONE
291C LAW STORAGE FOR DECAY AND HYSTERESIS
292C DO J=26,31
293C UVAR(I,J)=ZERO
294C ENDDO
295C ENDIF
296 ENDDO
297
298 IF(iflag==2)THEN
299 jacobi=.true.
300 iflag=0
301 ELSE
302 jacobi=.false.
303 ENDIF
304 IF (krecover==2) THEN
305 IF (time==zero) THEN
306 efac(1:nel) =zero
307 uvar(1:nel,32)=em20
308 ELSE
309 DO i=1,nel
310 uvar(i,32)=max(uvar(i,32),eint(i))
311 ! EFAC(I) =(ONE-HYSTER)*(ONE-(EINT(I)/UVAR(I,32))**BETA)
312 ! EFAC(I) =(ONE-HYSTER)*(ONE-PUI)
313 ! if eint = uvar --> Pui = 1 and efac = 0
314 ! if eint /= 0 --> Pui = exp(beta * ln(eint/uvar) )
315 ! if eint = 0 --> Pui = 0
316 IF(eint(i)/=uvar(i,32)) THEN
317 IF( (eint(i)/uvar(i,32))>zero ) THEN
318 pui = exp( beta*log( eint(i)/uvar(i,32) ) )
319 ELSE
320 pui = zero
321 ENDIF
322 ELSE
323 pui = one
324 ENDIF
325 efac(i) =(one-hyster)*( one- pui )
326 END DO
327 END IF
328 END IF
329C......................................................................
330
331C DEFINE STRETCH TENSOR
332
333
334 IF(iflag==0)THEN
335C PRINCIPAL STRAINS FORMULATION
336 mdir=3
337 DO i=1,nel
338 IF(total) THEN
339 av(i,1)=epsxx(i)
340 av(i,2)=epsyy(i)
341 av(i,3)=epszz(i)
342 av(i,4)=epsxy(i)*half
343 av(i,5)=epsyz(i)*half
344 av(i,6)=epszx(i)*half
345 ELSE
346 av(i,1)=depsxx(i)
347 av(i,2)=depsyy(i)
348 av(i,3)=depszz(i)
349 av(i,4)=depsxy(i)*half
350 av(i,5)=depsyz(i)*half
351 av(i,6)=depszx(i)*half
352 ENDIF
353 ENDDO
354 IF(jacobi) THEN
355 DO i=1,nel
356 a(1,1)=av(i,1)
357 a(2,2)=av(i,2)
358 a(3,3)=av(i,3)
359 a(2,1)=av(i,4)
360 a(3,2)=av(i,5)
361 a(3,1)=av(i,6)
362 CALL jacobiew(a,3,earj,dpra,nrot)
363 earv(i,1)=earj(1)
364 earv(i,2)=earj(2)
365 earv(i,3)=earj(3)
366 DO n=1,3
367 DO j=1,3
368 dirprv(i,n,j)=dpra(n,j)
369 ENDDO
370 ENDDO
371 ENDDO
372 ELSE
373C Eigenvalues needed to be calculated in double precision
374C for a simple precision executing
375 IF (iresp==1) THEN
376 CALL valpvecdp_v(av,earv,dirprv,nel)
377 ELSE
378 CALL valpvec_v(av,earv,dirprv,nel)
379 ENDIF
380 ENDIF
381C SET OLD PRINCIPAL DIRECTIONS
382 DO i=1,nel
383 m=0
384 DO n=1,3
385 DO l=1,3
386 m=m+1
387 dprao(l,n)=uvar(i,16+m)
388 ENDDO
389 ENDDO
390 m=0
391 DO n=1,3
392 DO j=1,3
393 m=m+1
394 dpra(n,j)=dirprv(i,n,j)
395 ENDDO
396 ENDDO
397
398C COMPARE NEW AND OLD PRINCIPAL DIRECTIONS AND TRANSFORM STORED
399C VALUES FROM OLD INTO NEW PRINCIPAL DIRECTIONS IF TOLERANCE EXCEEDED
400
401 CALL checkaxes(dprao,dpra,amin,amax,unchanged,tolerance)
402 IF(.NOT.unchanged) THEN
403 DO m=1,4
404 DO n=1,3
405 DO l=1,3
406 sn(n,l)=zero
407 IF(n==l) sn(n,l)= uvar(i,n+3*(m-1))
408 ENDDO
409 ENDDO
410C ROTATE STORED VALUES IN OLD PRINCIPAL DIRECTIONS INTO GLOBAL DIRECTIONS
411 ii = 1
412 jj = 1
413 ken = 1
414 CALL dreh(sn,dprao,ii,jj,ken)
415
416C ROTATE STORED VALUES IN GLOBAL DIRECTIONS INTO NEW PRINCIPAL DIRECTIONS
417 ii = 1
418 jj = 1
419 ken = 0
420 CALL dreh(sn,dpra ,ii,jj,ken)
421
422 DO n=1,3
423 uvar(i,n+3*(m-1)) = sn(n,n)
424 ENDDO
425 ENDDO
426C pour la partie decay et hysteresis
427 IF (beta>tiny.AND.krecover/=2) THEN
428 DO m=1,2
429 DO n=1,3
430 DO l=1,3
431 sn(n,l)=zero
432 IF(n==l) sn(n,l)= uvar(i,n+3*(m-1)+25)
433 ENDDO
434 ENDDO
435 ii = 1
436 jj = 1
437 ken = 1
438 CALL dreh(sn,dprao ,ii,jj,ken)
439
440 ii = 1
441 jj = 1
442 ken = 0
443 CALL dreh(sn,dpra ,ii,jj,ken)
444
445 DO n=1,3
446 uvar(i,n+3*(m-1)+25) = sn(n,n)
447 ENDDO
448 ENDDO
449 ENDIF
450 uvar(i,17)=dpra(1,1)
451 uvar(i,18)=dpra(2,1)
452 uvar(i,19)=dpra(3,1)
453 uvar(i,20)=dpra(1,2)
454 uvar(i,21)=dpra(2,2)
455 uvar(i,22)=dpra(3,2)
456 uvar(i,23)=dpra(1,3)
457 uvar(i,24)=dpra(2,3)
458 uvar(i,25)=dpra(3,3)
459 ENDIF
460
461 IF(total) THEN
462 ear(i,1)=earv(i,1)
463 ear(i,2)=earv(i,2)
464 ear(i,3)=earv(i,3)
465 ELSE
466 ecr(i,1)=earv(i,1)
467 ecr(i,2)=earv(i,2)
468 ecr(i,3)=earv(i,3)
469 ENDIF
470 ENDDO
471 ELSE
472
473C AVERAGE STRAIN FORMULATION
474
475 mdir=1
476 DO i=1,nel
477 IF(total) THEN
478 ear(i,1)= sign(one,epsxx(i)+epsyy(i)+epszz(i))*sqrt(
479 & (epsxx(i)-epsyy(i))**2+
480 & (epsyy(i)-epszz(i))**2+
481 & (epszz(i)-epsxx(i))**2+
482 & two*(epsxy(i)**2+epsyz(i)**2+epszx(i)**2))/sqrt(two)
483 ELSE
484 ecr(i,1)= sign(one,depsxx(i)+depsyy(i)+depszz(i))*sqrt(
485 & (depsxx(i)-depsyy(i))**2+
486 & (depsyy(i)-depszz(i))**2+
487 & (depszz(i)-depsxx(i))**2+
488 & ((depsxy(i))**2+(depsyz(i))**2+(depszx(i))**2))
489 & /sqrt(two)
490 ENDIF
491 ENDDO
492
493 ENDIF
494
495 IF(ismstr==10.OR.ismstr==12)THEN
496 DO i=1,nel
497 IF(off(i)==zero )THEN
498 ear(i,1)=zero
499 ear(i,2)=zero
500 ear(i,3)=zero
501 END IF
502 END DO
503 END IF
504C......................................................................
505 IF(timestep<=zero) THEN
506 dti=zero
507 ELSE
508 dti=one/timestep
509 ENDIF
510
511 DO j=1,mdir
512 DO i=1,nel
513 IF( total) THEN
514 IF ((ismstr==10.OR.ismstr==12).AND.ideac==0.AND.offg(i) <=one) THEN
515 ecr(i,j)=sqrt(ear(i,j)+one)-uvar(i,j)
516 ELSE
517 ecr(i,j)=ear(i,j)-uvar(i,j)
518 END IF
519 ELSE
520 ear(i,j)=ecr(i,j)+uvar(i,j)
521 END IF
522 ebr(i,j)=ecr(i,j)*dti
523 ENDDO
524C COMPUTE NOMINAL VALUES
525 DO i=1,nel
526C COMPUTE PRINCIPAL STRETCHES
527C LAMDA=E_NOMINAL+1=EXP(E_RADIOSS)
528 IF(ismstr==0.OR.ismstr==2.OR.ismstr==4) THEN
529 el(i,j)=exp(ear(i,j))
530C COMPUTE PRINCIPAL NOMINAL STRAIN(J) RATES FROM PRINCIPAL STRAIN(J) RATES
531C E_DOT_NOMINAL=EDOT_RADIOSS*EXP(E_RADIOSS)=EDOT_RADIOSS*LAMDA
532 ebn(i,j)=ebr(i,j)*el(i,j)
533 ecn(i,j)=el(i,j)*(one-exp(-ecr(i,j)))
534 ELSEIF(ismstr==10.OR.ismstr==12) THEN
535 IF (offg(i) <=one) THEN
536 el(i,j)=sqrt(ear(i,j)+one)
537 ELSE
538 el(i,j)=ear(i,j)+one
539 END IF
540 ebn(i,j)=ebr(i,j)
541 ecn(i,j)=ecr(i,j)
542 ELSE
543 el(i,j)=ear(i,j)+one
544 ebn(i,j)=ebr(i,j)
545 ecn(i,j)=ecr(i,j)
546 ENDIF
547C COMPUTE PRINCIPAL NOMINAL STRAINS FROM PRINCIPAL STRAINS
548 ean(i,j)=el(i,j)-one
549C STRAIN RATE FILTERING
550 ebn(i,j)=uvar(i,j+6)+(ebn(i,j)-uvar(i,j+6))*ratedamp
551 ENDDO
552 ENDDO
553
554C......................................................................
555
556C RELATIVE VOLUME
557 DO i=1,nel
558 pair(i)=zero
559 volumer(i)=el(i,1)*el(i,2)*el(i,3)
560 ENDDO
561
562C CONFINED AIR PRESSURE
563 DO i=1,nel
564 IF (kcompair==1.AND.volumer(i)<one) THEN
565 npcurve=ifunc(nfuncp)
566 IF (volumer(i)-phi>small) THEN
567 IF(npcurve/=0) THEN
568 pair(i)=-finter(npcurve,volumer(i),npf,tf,df)
569 pair(i)=pair(i)*pscale
570 ELSE
571 pair(i)=p0*(volumer(i)-one)/(volumer(i)-phi)
572 ENDIF
573 ELSE
574 pair(i)=uvar(i,16)
575 ENDIF
576C RELAXATION ON AIR PRESSURE
577 pair(i)=exp(-relaxp*time)*max(pair(i),-maxpres)
578 uvar(i,16)=pair(i)
579 ELSEIF (kcompair==2.AND.volumer(i)<one) THEN
580 npcurve=ifunc(nfuncp)
581 pair(i)=-finter(npcurve,volumer(i),npf,tf,df)
582 ENDIF
583 ENDDO
584
585C......................................................................
586
587C UPDATE TENSILE MODULE (INCREASES WITH COMPRESSION)
588
589C COMPUTE STRESS COMPONENTS VIA INTERPOLATION
590
591C LOOP ON PRINCIPAL DIRECTIONS
592 DO j=1,mdir
593 DO i=1,nel
594 psn(i,j)=zero
595C COMPRESSION IS POSITIVE!
596 strain(i,j)=-ean(i,j)
597 rate(i,j)=abs(ebn(i,j))
598 shift=zero
599 unloading=.false.
600 IF((ean(i,j)<zero.AND.ebn(i,j)>zero).OR.
601 & (ean(i,j)>zero.AND.ebn(i,j)<zero))unloading=.true.
602C COMPUTE ELASTIC MODULE FROM STATIC CURVE
603 psn1(i,j)=-alpha(1)*finter(ifunc(1),strain(i,j),npf,tf,df)
604C CURVE INTERPOLATION
605C EDOT0 = ACTUAL RATE
606C EDOTU = UNLOADING RATE
607C EDOTL = UPPER BOUND RATE
608C EDOTS = LOWERBOUND RATE
609 edot0(i,j)=rate(i,j)
610 IF(edot0(i,j)<edot(1)) edot0(i,j)=edot(1)
611 IF(unloading.AND.iunload/=0) THEN
612C CASE OF IUNLOAD/=0 (UNLOAD BETWEEN THE IUNLOAD CURVE (UPPER LIMIT)
613C AND THE FIRST CURVE (LOWER LIMIT)- STATIC- )
614 IF (abs(runload-edot(1))<em20) THEN
615 psn(i,j)=-funload*
616 & (finter(ifunc(nfuncul),strain(i,j),npf,tf,df1))
617 eyn(i,j)=funload*df1
618 visc(i,j)=zero
619 ELSE
620 edotu(i,j)=edot0(i,j)
621 edots(i,j)=edot(1)
622 edotl(i,j)=max(runload,edots(i,j))
623
624 IF(edotu(i,j)<edots(i,j))edotu(i,j)=edots(i,j)
625 IF(edotu(i,j)>edotl(i,j))edotu(i,j)=edotl(i,j)
626
627 psn1(i,j)=-alpha(1)*
628 & finter(ifunc(1),strain(i,j),npf,tf,df1)
629 eyn(i,j)=alpha(1)*df1
630 psn2(i,j)=-funload*
631 & (finter(ifunc(nfuncul),strain(i,j),npf,tf,df2))
632
633 ratio = (min(one,(edotu(i,j)-edots(i,j))/
634 & (max(tiny,edotl(i,j)-edots(i,j)))))
635
636 IF(ratio>zero) THEN
637 pui_1 = exp( exponas*log(ratio) )
638 ELSE
639 pui_1 = zero
640 ENDIF
641 IF(pui_1 < one) THEN
642 psn(i,j)=psn2(i,j)+(psn1(i,j)-psn2(i,j))*
643 & exp( exponbs*log( one - pui_1 ) )
644 ELSE
645 psn(i,j)=psn2(i,j)
646 ENDIF
647
648! PSN(I,J)=PSN2(I,J)+(PSN1(I,J)-PSN2(I,J))*
649! & EXP( EXPONBS*LOG( ONE-EXP( EXPONAS*LOG(RATIO) ) ) )
650
651! PSN(I,J)=PSN2(I,J)+(PSN1(I,J)-PSN2(I,J))*
652! & ( ONE-((MIN(ONE,(EDOTU(I,J)-EDOTS(I,J))/
653! & (MAX(TINY,EDOTL(I,J)-EDOTS(I,J)))))**EXPONAS))**EXPONBS
654
655! PSN(I,J)=PSN2(I,J)+(PSN1(I,J)-PSN2(I,J))*
656! & EXP( EXPONBS * LOG( ONE-EXP(EXPONAS*LOG(RATIO)) ) )
657! PSN(I,J)=PSN2(I,J)+(PSN1(I,J)-PSN2(I,J))*
658! & (ONE-(RATIO)**EXPONAS) ** EXPONBS
659
660
661 IF(abs(edotu(i,j)-edots(i,j))>tiny.AND.
662 & abs(psn(i,j)-psn1(i,j))>tiny)visc(i,j)=
663 & abs((psn(i,j)-psn1(i,j))/(edotu(i,j)-edots(i,j)))
664
665 visc(i,j)=min(visc(i,j),viscosity)
666 psn(i,j)=psn1(i,j)-visc(i,j)*(edotu(i,j)-edots(i,j))
667 ENDIF
668
669 ELSE
670 kun(i,j)=0
671 IF(unloading) kun(i,j)=nfunc1
672C IF(UNLOADING) SHIFT=TWO*(SIGN(1,STRAIN(I,J))-STRAIN(I,J))
673C CASE OF ONLY ONE CURVE INPUT I.E. NO STRAIN(J) RATE EFFECT LOADING
674C AND UNLOADING FOLLOW THE SAME PATH
675 IF (nfunc1==1) THEN
676 psn(i,j)=-alpha(1)*
677 & finter(ifunc(kun(i,j)+1),strain(i,j)+shift,npf,tf,df)
678 eyn(i,j)=alpha(1)*df
679 ELSE
680C CASE OF IUNLOAD==0 ( UNLOAD BETWEEN THE IUNLOAD CURVE (UPPER LIMIT)
681C AND THE FIRST CURVE (LOWER LIMIT )- STATIC-
682C IF UNLOADING CURVES ARE SPECIFIED THEY WILL BE USED FOR INTERPOLATION
683C OTHERWISE STARTER INITIALIZES THEM TO THE DEFAULT CURVE NUMBER I.E. THE
684C FIRST LOADING (STATIC) CURVE
685 DO l=1,nfunc1
686 IF(edot0(i,j)<=edot(l)) GOTO 10
687 ENDDO
688 l=nfunc1
689 edot0(i,j)=edot(l)
690 10 CONTINUE
691 k(i,j)=l
692 kf(i,j)=l+kun(i,j)
693
694 psn2(i,j)=-alpha(k(i,j))*
695 & finter(ifunc(kf(i,j)),strain(i,j)+shift,npf,tf,df2)
696 k1(i,j)=max(1,l-1)
697 kf1(i,j)=max(1,l-1)+kun(i,j)
698 edotl(i,j)=edot(k(i,j))
699 edots(i,j)=edot(k1(i,j))
700 psn1(i,j)=-alpha(k1(i,j))*
701 & finter(ifunc(kf1(i,j)),strain(i,j)+shift,npf,tf,df1)
702 IF(l==nfunc1) THEN
703 eyn(i,j)=alpha(l)*df2
704 ELSE
705 eyn(i,j)=alpha(k1(i,j))*df1
706 ENDIF
707
708 ratio = (min(one,(edot0(i,j)-edots(i,j))/
709 & (max(tiny,edotl(i,j)-edots(i,j)))))
710 IF(ratio>zero) THEN
711 pui_1 = exp( exponas*log(ratio) )
712 ELSE
713 pui_1 = zero
714 ENDIF
715 IF(pui_1 < one) THEN
716 psn(i,j)=psn2(i,j)+(psn1(i,j)-psn2(i,j))*
717 & exp(exponbs*log(one - pui_1 ) )
718 ELSE
719 psn(i,j)=psn2(i,j)
720 ENDIF
721! PSN(I,J)=PSN2(I,J)+(PSN1(I,J)-PSN2(I,J))*
722! & EXP(EXPONBS*LOG(ONE-EXP( EXPONAS*LOG(RATIO) ) ) )
723
724! PSN(I,J)=PSN2(I,J)+(PSN1(I,J)-PSN2(I,J))*
725! & (ONE-(RATIO)**EXPONAS) ** EXPONBS
726
727
728
729 IF(abs(edot0(i,j)-edots(i,j))>tiny.AND.
730 & abs(psn(i,j)-psn1(i,j))>tiny)visc(i,j)=
731 & abs((psn(i,j)-psn1(i,j))/(edot0(i,j)-edots(i,j)))
732
733 visc(i,j)=min(visc(i,j),viscosity)
734 psn(i,j)=psn1(i,j)-visc(i,j)*(edot0(i,j)-edots(i,j))
735 ENDIF
736 ENDIF
737
738 ENDDO
739
740C......................................................................
741
742C UNLOADING
743
744 DO i=1,nel
745 decay = uvar(i,j+28)
746C HYSTERESIS AND DECAY OCCURS ON COMPRESSION ONLY BUT REMAINS PERMANENT
747C IN TENSION
748 IF(ean(i,j)<zero) THEN
749C NO RECOVERY ON UNLOADING - ACCUMULATE STRAIN
750 IF(krecover==0.AND.ecn(i,j)<zero)
751 & uvar(i,j+25)=uvar(i,j+25)+abs(ecn(i,j))
752C RECOVERY ON UNLOADING - DOWN INCREMENT STRAIN
753 IF(krecover==1)uvar(i,j+25)=uvar(i,j+25)-ecn(i,j)
754 IF(ebn(i,j)<zero) THEN
755 decay = min(one,hyster*(one-exp(-beta*uvar(i,j+25))))
756 ENDIF
757 ENDIF
758C----loi70 Iflag=4 :
759 IF (krecover==2) decay = efac(i)
760C DECAY ON LOADING AND UNLOADING
761 IF(kdecay==0)
762 & psn(i,j)=psn(i,j)*(one-decay)
763C DECAY ON LOADING ONLY
764 IF(kdecay==1.AND.ecn(i,j)<zero)
765 & psn(i,j)=psn(i,j)*(one-decay)
766C DECAY ON UNLOADING ONLY
767 IF(kdecay==2.AND.ecn(i,j)>zero)
768 & psn(i,j)=psn(i,j)*(one-decay)
769 uvar(i,j+28)=decay
770C......................................................................
771 dsigma(i,j)=psn(i,j)-uvar(i,j+3)
772 IF(total) THEN
773C SECANT MODULE
774 tmp1 = eyn(i,j)
775 ELSE
776 IF(abs(ecn(i,j))>tiny)THEN
777 eyn(i,j)=abs(dsigma(i,j)/ecn(i,j))
778 ELSE
779 eyn(i,j)=uvar(i,j+9)
780 ENDIF
781 ENDIF
782 ENDDO
783
784 IF(.NOT.total) THEN
785 DO i=1,nel
786 IF(sign(one,(ebn(i,j)/(uvar(i,j+6)+tiny)))/=one)
787 & eyn(i,j)=uvar(i,j+9)
788 ENDDO
789
790 DO i=1,nel
791 eyn(i,j)=eyn(i,j)*theta+uvar(i,j+9)*(one-theta)
792 ENDDO
793 ENDIF
794
795
796 IF(itotal==0.OR.itotal==2) THEN
797C LINEAR IN TENSION
798 DO i=1,nel
799 IF(ean(i,j)>=zero) THEN
800 uvar(i,j+28)=zero
801 tmp1=exp(-lamda*(volumer(i)-1.+epsfin))
802 ei(i)=efinal+(e0-efinal)*(1-tmp1)
803 tmp2=lamda*(efinal-e0)*tmp1
804 eyn(i,j)=max(ei(i),tmp2)
805 IF(total) THEN
806 psn(i,j)=ei(i)*ean(i,j)
807 ELSE
808 psn(i,j)=uvar(i,j+3)+ei(i)*ecn(i,j)
809 ENDIF
810 uvar(i,j+12)=vt/ei(i)
811 ENDIF
812 ENDDO
813C----loi70 Iflag=4 :
814 IF (krecover==2) THEN
815 DO i=1,nel
816 IF(ean(i,j)>=zero) THEN
817 decay = efac(i)
818C DECAY ON LOADING AND UNLOADING
819 IF(kdecay==0) psn(i,j)=psn(i,j)*(one-decay)
820C DECAY ON LOADING ONLY
821 IF(kdecay==1.AND.ecn(i,j)<zero)
822 & psn(i,j)=psn(i,j)*(one-decay)
823C DECAY ON UNLOADING ONLY
824 IF(kdecay==2.AND.ecn(i,j)>zero)
825 & psn(i,j)=psn(i,j)*(one-decay)
826 END if!(EAN(I,J)>=ZERO) THEN
827 END DO
828 END IF !(KRECOVER==2) THEN
829 ELSEIF(itotal==-1) THEN
830C MODULE MAXMAL IN TENSION ONE CYCLE LATER (TOTAL)
831 DO i=1,nel
832 IF(ean(i,j)>0.) THEN
833 uvar(i,j+28)=zero
834 ei(i)=max(uvar(i,10),uvar(i,11),uvar(i,12))
835 eyn(i,j)=ei(i)
836 psn(i,j)=ei(i)*ean(i,j)
837 uvar(i,j+12)=vt/ei(i)
838 ENDIF
839 ENDDO
840 ELSEIF(itotal==-2) THEN
841C MODULE MAXMAL IN TENSION ONE CYCLE LATER (INCR)
842 DO i=1,nel
843 IF(ean(i,j)>zero) THEN
844 uvar(i,j+28)=zero
845 ei(i)=max(uvar(i,10),uvar(i,11),uvar(i,12))
846 eyn(i,j)=ei(i)
847 psn(i,j)=uvar(i,j+3)+ei(i)*ecn(i,j)
848 uvar(i,j+12)=vt/ei(i)
849 ENDIF
850 ENDDO
851 ENDIF
852
853C END OF LOOP ON MDIR
854 ENDDO
855C BEGIN OF LOOP ON MDIR
856 IF (imsta>=1) THEN
857 DO i=1,nel
858 sigmax(i)=-(min(psn(i,1),psn(i,2),psn(i,3))-
859 . max(psn(i,1),psn(i,2),psn(i,3)))
860 ENDDO
861 DO j=1,mdir
862 DO i=1,nel
863 esecant(i,j)=0.4*abs(psn(i,j))/max(tiny,abs(strain(i,j)))
864 IF (esecant(i,j)<=sigmax(i)) THEN
865 tmp1=0.2*(sigmax(i)-esecant(i,j))
866 psn(i,j)=psn(i,j)+tmp1*ean(i,j)
867 eyn(i,j)=max(eyn(i,j),(one+tmp1)*esecant(i,j))
868 ENDIF
869 ENDDO
870 ENDDO
871 ENDIF
872C
873 IF ((ismstr==10.OR.ismstr==12).AND.ideac==0) THEN
874 DO j=1,mdir
875 DO i=1,nel
876 tmp0=max(eyn(i,j),uvar(i,j+9))
877 IF (offg(i) <=one) THEN
878 uvar(i,j )=sqrt(ear(i,j)+one)
879 ELSE
880 uvar(i,j )=ear(i,j)
881 END IF
882 uvar(i,j+3)=psn(i,j)
883 uvar(i,j+6)=ebn(i,j)
884 uvar(i,j+9)=eyn(i,j)
885 eyn(i,j)=tmp0
886 ENDDO
887C END OF LOOP ON MDIR
888 ENDDO
889 ELSE
890 DO j=1,mdir
891 DO i=1,nel
892 tmp0=max(eyn(i,j),uvar(i,j+9))
893 uvar(i,j )=ear(i,j)
894 uvar(i,j+3)=psn(i,j)
895 uvar(i,j+6)=ebn(i,j)
896 uvar(i,j+9)=eyn(i,j)
897 eyn(i,j)=tmp0
898 ENDDO
899C END OF LOOP ON MDIR
900 ENDDO
901 END IF !(ISMSTR==10.AND.IDEAC==0) THEN
902C---------for output EPSD :
903 DO i=1,nel
904 epsd(i)=max(abs(ebn(i,1)),abs(ebn(i,2)),abs(ebn(i,3)))
905 ENDDO
906
907 IF(iflag==1) THEN
908C OCTAHEDRAL STRAIN (AVERAGE) FORMULATION
909 DO i=1,nel
910 emax(i)=ei(i)
911 d11(i)=emax(i)*(one-vt)/(one+vt)/(one-two*vt)
912 d12(i)=emax(i)*vt/(one+vt)/(one-two*vt)
913 d44(i)=emax(i)/two/(one+vt)
914 ENDDO
915 DO i=1,nel
916 IF(total) THEN
917 signxx(i)=d11(i)*epsxx(i)+d12(i)*epsyy(i)+d12(i)*epszz(i)
918 signyy(i)=d12(i)*epsxx(i)+d11(i)*epsyy(i)+d12(i)*epszz(i)
919 signzz(i)=d12(i)*epsxx(i)+d12(i)*epsyy(i)+d11(i)*epszz(i)
920 signxy(i)=d44(i)*epsxy(i)
921 signyz(i)=d44(i)*epsyz(i)
922 signzx(i)=d44(i)*epszx(i)
923 soundsp(i) = sqrt(d11(i)/rho0(i))
924 ELSE
925 signxx(i)=sigoxx(i)+
926 & d11(i)*depsxx(i)+d12(i)*depsyy(i)+d12(i)*depszz(i)
927 signyy(i)=sigoyy(i)+
928 & d12(i)*depsxx(i)+d11(i)*depsyy(i)+d12(i)*depszz(i)
929 signzz(i)=sigozz(i)+
930 & d12(i)*depsxx(i)+d12(i)*depsyy(i)+d11(i)*depszz(i)
931 signxy(i)=sigoxy(i)+d44(i)*depsxy(i)
932 signyz(i)=sigoyz(i)+d44(i)*depsyz(i)
933 signzx(i)=sigozx(i)+d44(i)*depszx(i)
934 soundsp(i) = sqrt(d11(i)/rho0(i))
935 ENDIF
936 viscmax(i)=visc(i,1)
937
938 ENDDO
939 RETURN
940 ENDIF
941
942 DO i=1,nel
943 IF(vt+vc<=two*tiny) THEN
944 IF(total) THEN
945 psc(i,1)=psn(i,1)
946 psc(i,2)=psn(i,2)
947 psc(i,3)=psn(i,3)
948 ELSE
949 psc(i,1)=eyn(i,1)*ecn(i,1)
950 psc(i,2)=eyn(i,2)*ecn(i,2)
951 psc(i,3)=eyn(i,3)*ecn(i,3)
952 ENDIF
953 ELSE
954 e12(i)=(ean(i,1)+ean(i,2))/two
955 e23(i)=(ean(i,2)+ean(i,3))/two
956 e31(i)=(ean(i,3)+ean(i,1))/two
957 v12(i)=vc+(vt-vc)*(one-exp(-rv*abs(e12(i))))*
958 & (sign(one,e12(i))+one)/two
959 v23(i)=vc+(vt-vc)*(one-exp(-rv*abs(e23(i))))*
960 & (sign(one,e23(i))+one)/two
961 v31(i)=vc+(vt-vc)*(one-exp(-rv*abs(e31(i))))*
962 & (sign(one,e31(i))+one)/two
963 IF(total) THEN
964 detc(i)=one-v23(i)*v23(i)-v31(i)*v31(i)-v12(i)*v12(i)-
965 & two*v12(i)*v31(i)*v23(i)
966 a11(i)=(one-v23(i)*v23(i))
967 a12(i)=(v12(i)+v23(i)*v31(i))
968 a13(i)=(v31(i)+v23(i)*v12(i))
969 a22(i)=(one-v31(i)*v31(i))
970 a23(i)=(v23(i)+v31(i)*v12(i))
971 a33(i)=(one-v12(i)*v12(i))
972 psc(i,1)=(a11(i)*psn(i,1)+a12(i)*psn(i,2)+a13(i)*psn(i,3))/
973 & detc(i)
974 psc(i,2)=(a12(i)*psn(i,1)+a22(i)*psn(i,2)+a23(i)*psn(i,3))/
975 & detc(i)
976 psc(i,3)=(a13(i)*psn(i,1)+a23(i)*psn(i,2)+a33(i)*psn(i,3))/
977 & detc(i)
978 ELSE
979 uvar(i,13)=theta*v23(i)/ei(i)+(one-theta)*uvar(i,13)
980 uvar(i,14)=theta*v31(i)/ei(i)+(one-theta)*uvar(i,14)
981 uvar(i,15)=theta*v12(i)/ei(i)+(one-theta)*uvar(i,15)
982 detc(i)=one/(eyn(i,1)*eyn(i,2)*eyn(i,3))-
983 & uvar(i,13)*uvar(i,13)/eyn(i,1)-
984 & uvar(i,14)*uvar(i,14)/eyn(i,2)-
985 & uvar(i,15)*uvar(i,15)/eyn(i,3)-
986 & 2*uvar(i,13)*uvar(i,14)*uvar(i,15)
987 a11(i)=one/(eyn(i,2)*eyn(i,3))-uvar(i,13)*uvar(i,13)
988 a12(i)=uvar(i,15)/eyn(i,3)+uvar(i,13)*uvar(i,14)
989 a13(i)=uvar(i,14)/eyn(i,2)+uvar(i,13)*uvar(i,15)
990 a22(i)=one/(eyn(i,1)*eyn(i,3))-uvar(i,14)*uvar(i,14)
991 a23(i)=uvar(i,13)/eyn(i,1)+uvar(i,14)*uvar(i,15)
992 a33(i)=one/(eyn(i,1)*eyn(i,2))-uvar(i,15)*uvar(i,15)
993C COMPUTE STRESS INCREMENT IN PRINCIPAL DIRECTIONS
994 psc(i,1)=((a11(i)*ecn(i,1)+a12(i)*ecn(i,2)+a13(i)*ecn(i,3))
995 & /detc(i))
996 psc(i,2)=((a12(i)*ecn(i,1)+a22(i)*ecn(i,2)+a23(i)*ecn(i,3))
997 & /detc(i))
998 psc(i,3)=((a13(i)*ecn(i,1)+a23(i)*ecn(i,2)+a33(i)*ecn(i,3))
999 & /detc(i))
1000 ENDIF
1001 ENDIF
1002
1003 DO j=1,3
1004 IF(off(i)==zero.OR.psn(i,j)>abs(tensioncut))THEN
1005 psc(i,1)=zero
1006 psc(i,2)=zero
1007 psc(i,3)=zero
1008 off(i)=zero
1009 ENDIF
1010 ENDDO
1011C COMPUTE CAUCHY STRESS INCREMENT IN PRINCIPAL DIRECTIONS
1012CSIGMA_CAUCHY(I) = LAMDA(I) * SIGMA_NOMINAL(I) / RELATIVE VOLUME
1013 IF(ismstr==0.OR.ismstr==2.OR.ismstr==4) THEN
1014C CAUCHY MODULE FOR TIME STEP
1015 den1=el(i,2)*el(i,3)
1016 den2=el(i,3)*el(i,1)
1017 den3=el(i,1)*el(i,2)
1018 psc(i,1)=psc(i,1)/den1
1019 psc(i,2)=psc(i,2)/den2
1020 psc(i,3)=psc(i,3)/den3
1021 eyn(i,1)=eyn(i,1)/den1
1022 eyn(i,2)=eyn(i,2)/den2
1023 eyn(i,3)=eyn(i,3)/den3
1024 ELSEIF(ismstr==10.OR.(ismstr==12.AND.offg(i)<=one)) THEN
1025C------------directly to Chauchy-------
1026 den1=el(i,1)/volumer(i)
1027 den2=el(i,2)/volumer(i)
1028 den3=el(i,3)/volumer(i)
1029 psc(i,1)=psc(i,1)*den1
1030 psc(i,2)=psc(i,2)*den2
1031 psc(i,3)=psc(i,3)*den3
1032 eyn(i,1)=eyn(i,1)*den1
1033 eyn(i,2)=eyn(i,2)*den2
1034 eyn(i,3)=eyn(i,3)*den3
1035 ENDIF
1036 ENDDO
1037 IF (kcompair==2) THEN
1038 DO i=1,nel
1039 tmp0=volumer(i)
1040 tmp3=min(el(i,1),el(i,2),el(i,3))
1041 IF (tmp0<one.AND.tmp3<one
1042 & .AND.tmp3>tmp0.AND.abs(tmp0-tmp3)>em6) THEN
1043 IF(volumer(i)>zero) THEN
1044 tmp2=exp((1./3.)*log(volumer(i)))-volumer(i)
1045 ELSE
1046 tmp2=zero
1047 ENDIF
1048 tmp4=(tmp3-tmp0)/tmp2
1049 aa(i)=min(one,tmp4)
1050 ELSE
1051 aa(i)=zero
1052 ENDIF
1053 ENDDO
1054 DO j=1,3
1055 DO i=1,nel
1056 sigprv(i,j)=psc(i,j)+aa(i)*(pair(i)-psc(i,j))
1057 ENDDO
1058 ENDDO
1059 ELSE
1060 DO j=1,3
1061 DO i=1,nel
1062 sigprv(i,j)=psc(i,j)+pair(i)
1063 ENDDO
1064 ENDDO
1065 ENDIF
1066C TRANSFORM FROM PRINCIPAL TO GLOBAL DIRECTIONS
1067 DO i=1,nel
1068 signxx(i) = dirprv(i,1,1)*dirprv(i,1,1)*sigprv(i,1)
1069 & + dirprv(i,1,2)*dirprv(i,1,2)*sigprv(i,2)
1070 & + dirprv(i,1,3)*dirprv(i,1,3)*sigprv(i,3)
1071 signyy(i) = dirprv(i,2,2)*dirprv(i,2,2)*sigprv(i,2)
1072 & + dirprv(i,2,3)*dirprv(i,2,3)*sigprv(i,3)
1073 & + dirprv(i,2,1)*dirprv(i,2,1)*sigprv(i,1)
1074 signzz(i) = dirprv(i,3,3)*dirprv(i,3,3)*sigprv(i,3)
1075 & + dirprv(i,3,1)*dirprv(i,3,1)*sigprv(i,1)
1076 & + dirprv(i,3,2)*dirprv(i,3,2)*sigprv(i,2)
1077 signxy(i) = dirprv(i,1,1)*dirprv(i,2,1)*sigprv(i,1)
1078 & + dirprv(i,1,2)*dirprv(i,2,2)*sigprv(i,2)
1079 & + dirprv(i,1,3)*dirprv(i,2,3)*sigprv(i,3)
1080 signyz(i) = dirprv(i,2,2)*dirprv(i,3,2)*sigprv(i,2)
1081 & + dirprv(i,2,3)*dirprv(i,3,3)*sigprv(i,3)
1082 & + dirprv(i,2,1)*dirprv(i,3,1)*sigprv(i,1)
1083 signzx(i) = dirprv(i,3,3)*dirprv(i,1,3)*sigprv(i,3)
1084 & + dirprv(i,3,1)*dirprv(i,1,1)*sigprv(i,1)
1085 & + dirprv(i,3,2)*dirprv(i,1,2)*sigprv(i,2)
1086 ENDDO
1087
1088 IF(.NOT.total) THEN
1089 DO i=1,nel
1090C ADD INCREMENT OF STRESS
1091 signxx(i)=signxx(i)+sigoxx(i)
1092 signyy(i)=signyy(i)+sigoyy(i)
1093 signzz(i)=signzz(i)+sigozz(i)
1094 signxy(i)=signxy(i)+sigoxy(i)
1095 signyz(i)=signyz(i)+sigoyz(i)
1096 signzx(i)=signzx(i)+sigozx(i)
1097 ENDDO
1098 ENDIF
1099C SOUNDSPEED
1100 IF(iflag==0) THEN
1101 DO i=1,nel
1102 emax(i)=max(ei(i),eyn(i,1),eyn(i,2),eyn(i,3))
1103 ENDDO
1104 ENDIF
1105
1106 IF (imsta==2) THEN
1107 DO i=1,nel
1108 epsxy(i) = dirprv(i,1,1)*dirprv(i,2,1)*ean(i,1)
1109 & + dirprv(i,1,2)*dirprv(i,2,2)*ean(i,2)
1110 & + dirprv(i,1,3)*dirprv(i,2,3)*ean(i,3)
1111 epsyz(i) = dirprv(i,2,2)*dirprv(i,3,2)*ean(i,2)
1112 & + dirprv(i,2,3)*dirprv(i,3,3)*ean(i,3)
1113 & + dirprv(i,2,1)*dirprv(i,3,1)*ean(i,1)
1114 epszx(i) = dirprv(i,3,3)*dirprv(i,1,3)*ean(i,3)
1115 & + dirprv(i,3,1)*dirprv(i,1,1)*ean(i,1)
1116 & + dirprv(i,3,2)*dirprv(i,1,2)*ean(i,2)
1117 ENDDO
1118 DO i=1,nel
1119 esecant(i,1)=0.5*abs(signxy(i))/max(tiny,abs(epsxy(i)))
1120 esecant(i,2)=0.5*abs(signyz(i))/max(tiny,abs(epsyz(i)))
1121 esecant(i,3)=0.5*abs(signzx(i))/max(tiny,abs(epszx(i)))
1122 sigmax(i)=max(0.5*ei(i),sigmax(i))
1123 IF (esecant(i,1)<=sigmax(i)) THEN
1124 tmp1=0.1*(sigmax(i)-esecant(i,1))
1125 signxy(i)=signxy(i)+tmp1*epsxy(i)
1126 ENDIF
1127 IF (esecant(i,2)<=sigmax(i)) THEN
1128 tmp1=0.1*(sigmax(i)-esecant(i,2))
1129 signyz(i)=signyz(i)+tmp1*epsyz(i)
1130 ENDIF
1131 IF (esecant(i,3)<=sigmax(i)) THEN
1132 tmp1=0.1*(sigmax(i)-esecant(i,3))
1133 signzx(i)=signzx(i)+tmp1*epszx(i)
1134 ENDIF
1135 ENDDO
1136 ENDIF
1137 DO i=1,nel
1138 kkk=emax(i)/(three*(one-two*max(vc,vt)))
1139 ggg=emax(i)/(two*(one+max(vc,vt)))
1140 soundsp(i)=sqrt((kkk+four_over_3*ggg)/rho0(i))
1141 viscmax(i)=max(visc(i,1),visc(i,2),visc(i,3))
1142 ENDDO
1143C
1144 IF (impl_s > 0 .OR. ihet > 1) THEN
1145C-----choix entre EYN(I,j) et ESECANT(I,j) ----
1146 IF(iflag/=0) THEN
1147 DO i=1,nel
1148 emax(i)=max(ei(i),eyn(i,1),eyn(i,2),eyn(i,3))
1149 ENDDO
1150 ENDIF
1151 DO i=1,nel
1152 eymin=min(eyn(i,1),eyn(i,2),eyn(i,3))
1153 eyav = third*(eyn(i,1)+eyn(i,2)+eyn(i,3))
1154 et(i)= eymin/emax(i)
1155 ENDDO
1156C---overestimated ET but stable for not very high E0
1157 IF (time==zero) THEN
1158 tmp1=-alpha(1)*finter(ifunc(1),em20,npf,tf,df)
1159 tmp2 = alpha(1)*df/e0
1160 IF (tmp2>0.05) THEN
1161 uvar(1:nel,33) = -one
1162 ELSE
1163 uvar(1:nel,33) =et(1:nel)
1164 END IF
1165 ELSEIF (nel>0.AND.uvar(1,33)<zero) THEN
1166 et(1:nel) = one
1167 ELSE
1168 tmp1 = 0.75
1169 et(1:nel)=tmp1*et(1:nel)+(one-tmp1)*uvar(1:nel,33)
1170 uvar(1:nel,33) =et(1:nel)
1171 END IF !(TIME==ZERO) THEN
1172 END IF !(IMPL_S > 0 .OR. IHET > 1) THEN
1173C
1174 RETURN
1175
1176 END
1177
1178!||====================================================================
1179!|| checkaxes ../engine/source/materials/mat/mat038/sigeps38.F
1180!||--- called by ------------------------------------------------------
1181!|| sigeps38 ../engine/source/materials/mat/mat038/sigeps38.F
1182!||====================================================================
1183 SUBROUTINE checkaxes(TRAN1,TRAN2,AMIN,AMAX,UNCHANGED,TOL)
1184C----------------------------------------------------------------
1185 IMPLICIT NONE
1186#include "constant.inc"
1187
1188 my_real
1189 . tran1(3,3),tran2(3,3),ang(3),amax,tol,amin
1190 INTEGER I,J,K
1191 LOGICAL UNCHANGED
1192C----------------------------------------------------------------
1193 AMAX=zero
1194 unchanged=.false.
1195 DO j=1,3
1196 ang(j)=zero
1197 DO k=1,3
1198 ang(j)=ang(j)+tran1(k,j)*tran2(k,j)
1199 ENDDO
1200 ang(j)=abs(abs(ang(j))-one)
1201 IF (ang(j)>amax) amax=ang(j)
1202 ENDDO
1203 IF (amax<tol) unchanged=.true.
1204 amin=ang(1)
1205 RETURN
1206 END
my_real
#define my_real
Definition cppsort.cpp:32
if
if(complex_arithmetic) id
sigeps38
subroutine sigeps38(nel, nuparam, nuvar, nfunc, ifunc, npf, tf, time, timestep, uparam, rho0, rho, volume, eint, epspxx, epspyy, epspzz, epspxy, epspyz, epspzx, depsxx, depsyy, depszz, depsxy, depsyz, depszx, epsxx, epsyy, epszz, epsxy, epsyz, epszx, sigoxx, sigoyy, sigozz, sigoxy, sigoyz, sigozx, signxx, signyy, signzz, signxy, signyz, signzx, sigvxx, sigvyy, sigvzz, sigvxy, sigvyz, sigvzx, soundsp, viscmax, uvar, off, ismstr, mfxx, mfxy, mfxz, mfyx, mfyy, mfyz, mfzx, mfzy, mfzz, et, ihet, offg, epsd)
Definition sigeps38.F:50
checkaxes
subroutine checkaxes(tran1, tran2, amin, amax, unchanged, tol)
Definition sigeps38.F:1184
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
dreh
subroutine dreh(b, tr, ii, jj, ken)
Definition matrix.F:161
jacobiew
subroutine jacobiew(a, n, ew, ev, nrot)
Definition matrix.F:29
valpvecdp_v
subroutine valpvecdp_v(sig, val, vec, nel)
Definition sigeps33.F:2231
valpvec_v
subroutine valpvec_v(sig, val, vec, nel)
Definition sigeps33.F:1902