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