OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps38.F File Reference
#include "implicit_f.inc"
#include "mvsiz_p.inc"
#include "scr05_c.inc"

Go to the source code of this file.

Functions/Subroutines

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, uvarbid, ismstr, mfxx, mfxy, mfxz, mfyx, mfyy, mfyz, mfzx, mfzy, mfzz)
subroutine checkaxes (tran1, tran2, amin, amax, unchanged, tol)

Function/Subroutine Documentation

◆ checkaxes()

subroutine checkaxes ( tran1,
tran2,
amin,
amax,
logical unchanged,
tol )

Definition at line 987 of file sigeps38.F.

988C----------------------------------------------------------------
989 IMPLICIT NONE
990#include "constant.inc"
991 my_real
992 . tran1(3,3),tran2(3,3),ang(3),amax,tol,amin
993 INTEGER I,J,K
994 LOGICAL UNCHANGED
995C----------------------------------------------------------------
996 amax=zero
997 unchanged=.false.
998 DO j=1,3
999 ang(j)=zero
1000 DO k=1,3
1001 ang(j)=ang(j)+tran1(k,j)*tran2(k,j)
1002 ENDDO
1003 ang(j)=abs(abs(ang(j))-one)
1004 IF (ang(j)>amax) amax=ang(j)
1005 ENDDO
1006 IF (amax<tol) unchanged=.true.
1007 amin=ang(1)
1008 RETURN
my_real
#define my_real
Definition cppsort.cpp:32

◆ sigeps38()

subroutine sigeps38 ( integer nel,
integer nuparam,
integer nuvar,
integer nfunc,
integer, dimension(nfunc) ifunc,
integer, dimension(*) 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,
uvarbid,
integer ismstr,
mfxx,
mfxy,
mfxz,
mfyx,
mfyy,
mfyz,
mfzx,
mfzy,
mfzz )

Definition at line 35 of file sigeps38.F.

48C-----------------------------------------------
49C I M P L I C I T T Y P E S
50C-----------------------------------------------
51#include "implicit_f.inc"
52C-----------------------------------------------
53C G l o b a l P a r a m e t e r s
54C-----------------------------------------------
55#include "mvsiz_p.inc"
56#include "scr05_c.inc"
57C----------------------------------------------------------------
58C I N P U T A R G U M E N T S
59C----------------------------------------------------------------
60 INTEGER NEL, NUPARAM, NUVAR,ISMSTR
61 my_real
62 . uvarbid
63 my_real
64 . time , timestep , uparam(nuparam),
65 . rho(nel), rho0(nel), volume(nel), eint(nel),
66 . epspxx(nel), epspyy(nel), epspzz(nel),
67 . epspxy(nel), epspyz(nel), epspzx(nel),
68 . depsxx(nel), depsyy(nel), depszz(nel),
69 . depsxy(nel), depsyz(nel), depszx(nel),
70 . epsxx(nel), epsyy(nel), epszz(nel),
71 . epsxy(nel), epsyz(nel), epszx(nel),
72 . sigoxx(nel), sigoyy(nel), sigozz(nel),
73 . sigoxy(nel), sigoyz(nel), sigozx(nel),
74 . mfxx(nel) , mfxy(nel), mfxz(nel),
75 . mfyx(nel) , mfyy(nel), mfyz(nel),
76 . mfzx(nel) , mfzy(nel), mfzz(nel)
77C----------------------------------------------------------------
78C O U T P U T A R G U M E N T S
79C----------------------------------------------------------------
80 my_real
81 . signxx(nel), signyy(nel), signzz(nel),
82 . signxy(nel), signyz(nel), signzx(nel),
83 . sigvxx(nel), sigvyy(nel), sigvzz(nel),
84 . sigvxy(nel), sigvyz(nel), sigvzx(nel),
85 . soundsp(nel), viscmax(nel),
86 . den1,den2,den3,den1d,den2d,
87 . den3d,old1,old2,old3,old4,old5,old6
88C----------------------------------------------------------------
89C I N P U T O U T P U T A R G U M E N T S
90C----------------------------------------------------------------
91 my_real
92 . uvar(nel,nuvar), off(nel)
93C----------------------------------------------------------------
94C VARIABLES FOR FUNCTION INTERPOLATION
95C----------------------------------------------------------------
96 INTEGER NPF(*), NFUNC, IFUNC(NFUNC)
97 my_real
98 . finter,tf(*)
99 EXTERNAL finter
100C----------------------------------------------------------------
101C L O C A L V A R I B L E S
102C----------------------------------------------------------------
103 INTEGER MFUNC,IUNLOAD
104 INTEGER NUPARAM0
105 INTEGER I,J,L,N,M,MDIR,NROT
106 INTEGER K(MVSIZ,3),K1(MVSIZ,3),KF(MVSIZ,3)
107 INTEGER KF1(MVSIZ,3),KUN(MVSIZ,3)
108 INTEGER IFLAG,ITOTAL,IMSTA
109 INTEGER NFUNC1,NFUNCUL,NFUNCP,NPCURVE,KCOMPAIR
110 INTEGER KRECOVER,KDECAY
111 LOGICAL TOTAL,UNCHANGED,JACOBI,UNLOADING
112 INTEGER IND1,LD
113 INTEGER II,JJ,KEN
114C REAL
115 my_real
116 . strain(mvsiz,3),rate(mvsiz,3),ratem(mvsiz),
117 . visc(mvsiz,3),
118 . dsigma(mvsiz,3),decay,tensioncut,
119 . amin,amax,tolerance,lamda,efinal,epsfin,efactor,
120 . a(3,3),sn(3,3),earj(3),
121 . psc(mvsiz,3),ear(mvsiz,3),
122 . ebr(mvsiz,3),ecr(mvsiz,3),
123 . psn(mvsiz,3),ean(mvsiz,3),
124 . ebn(mvsiz,3),ecn(mvsiz,3),
125 . el(mvsiz,3),
126 . dpra(3,3),dprao(3,3),dprad(3,3),
127 . av(6,mvsiz),earv(3,mvsiz),dirprv(3,3,mvsiz),
128 . sigprv(3,mvsiz),ei(mvsiz),
129 . e0,vt,vc,rv,kkk,ggg,
130 . beta,hyster,ratedamp,theta,
131 . p0,relaxp,maxpres,phi,gamma,
132 . pair(mvsiz),volumer(mvsiz),
133 . df,df1,df2,
134 . viscosity,
135 . alpha(5),edot(5),
136 . psn1(mvsiz,3),psn2(mvsiz,3),
137 . edot0(mvsiz,3),edots(mvsiz,3),
138 . edotl(mvsiz,3),edotu(mvsiz,3),
139 . exponas,exponbs,funload,runload,
140 . e12(mvsiz),e23(mvsiz),e31(mvsiz),
141 . a11(mvsiz),a12(mvsiz),a13(mvsiz),a22(mvsiz),
142 . a23(mvsiz),a33(mvsiz),detc(mvsiz),
143 . v12(mvsiz),v23(mvsiz),v31(mvsiz),
144 . d11(mvsiz),d12(mvsiz),d44(mvsiz),
145 . emax(mvsiz),eyn(mvsiz,3),
146 . huge,small,tiny,shift,pscale
147c . ESECANT(MVSIZ,3),ETANGENT(MVSIZ,3),EELASTIC(MVSIZ,3),
148 my_real
149 . ar1(6),dar1(6),ar2(6),dar2(6),dar3(3),earjd(3),
150 . tmp0,tmp1,tmp2,tmp3,tmp4,tmp5,tmp6,ad(3,3),
151 . aa(mvsiz),sigmax(mvsiz),esecant(mvsiz,3),
152 . tmp0p,tmp1p,tmp2p,tmp3p,tmp4p,tmp5p,tmp6p, dti
153C=======================================================================
154 huge = ep30
155 small= em3
156 tiny = em30
157C......................................................................
158C SET INITIAL MATERIAL CONSTANTS
159
160 nuparam0= uparam(1)
161
162 e0 = uparam(2)
163 vt = uparam(3)
164 vc = uparam(4)
165 rv = uparam(5)
166 iflag = uparam(6)
167 itotal = uparam(7)
168
169 beta = uparam(8)
170 hyster = uparam(9)
171 ratedamp = uparam(10)
172 krecover = uparam(11)
173 kdecay = uparam(12)
174 theta = uparam(13)
175c THETA = RATEDAMP
176
177 kcompair = uparam(14)
178 p0 = uparam(15)
179 gamma = uparam(16)
180 relaxp = uparam(17)
181 maxpres = uparam(18)
182 phi = uparam(19)
183
184 iunload = uparam(20)
185 funload = uparam(21)
186 runload = uparam(22)
187 exponas = uparam(23)
188 exponbs = uparam(24)
189
190 mfunc = uparam(25)
191 imsta = uparam(26)
192 tensioncut= uparam(27)
193
194 efinal = uparam(28)
195 epsfin = uparam(29)
196 lamda = uparam(30)
197 viscosity = uparam(31)
198 tolerance = uparam(32)
199 pscale = uparam(33)
200 nfunc1=(nfunc-2)/2
201* unloading function number
202 nfuncul=nfunc-1
203* function number for enclosed air pressure
204 nfuncp=nfunc
205
206* set strain rates and scale factors for curves
207 DO i=1,nfunc1
208 edot(i) =uparam(nuparam-nfunc1+i)
209 alpha(i)=uparam(nuparam-nfunc1*2+i)
210 ENDDO
211
212* initialize e-module
213 total=.true.
214C ITOTAL=0 ; TOTAL - LINEAR IN TENSION
215C ITOTAL=1 ; TOTAL - NONLINEAR IN TENSION
216C ITOTAL=2 ; INCREMENTAL - LINEAR IN TENSION
217C ITOTAL=3 ; INCREMENTAL - NONLINEAR IN TENSION
218 IF(itotal==2.OR.itotal==3) total=.false.
219C +++ ([F]=[M_F]+[1]) -----
220 IF(ismstr>=10) THEN
221C--------- [B]=[F][F]^t strain-----
222c DO I=1,NEL
223c EPSXX(I)=MFXX(I)*(TWO+MFXX(I))+
224c . MFXY(I)*MFXY(I)+MFXZ(I)*MFXZ(I)
225c EPSYY(I)=MFYY(I)*(TWO+MFYY(I))+
226c . MFYX(I)*MFYX(I)+MFYZ(I)*MFYZ(I)
227c EPSZZ(I)=MFZZ(I)*(TWO+MFZZ(I))+
228c . MFZX(I)*MFZX(I)+MFZY(I)*MFZY(I)
229c EPSXY(I)=TWO*(MFXY(I)+MFYX(I)+MFXX(I)*MFYX(I)+
230c . MFXY(I)*MFYY(I)+MFXZ(I)*MFYZ(I))
231c EPSZX(I)=TWO*(MFXZ(I)+MFZX(I)+MFXX(I)*MFZX(I)+
232c . MFXY(I)*MFZY(I)+MFXZ(I)*MFZZ(I))
233c EPSYZ(I)=TWO*(MFZY(I)+MFYZ(I)+MFZX(I)*MFYX(I)+
234c . MFZY(I)*MFYY(I)+MFZZ(I)*MFYZ(I))
235c ENDDO
236 total=.true.
237 ENDIF
238
239C......................................................................
240C INITIALIZE
241 DO i=1,nel
242 ei(i)=e0
243 DO j=1,3
244 visc(i,j)=zero
245 ENDDO
246C IF(TIME==ZERO) THEN
247C EYN(I,1) = E0
248C EYN(I,2) = E0
249C EYN(I,3) = E0
250C STRAIN(1-3),STRESS(4-6),STRAIN RATE(7-9)
251C DO J=1,9
252C UVAR(I,J)=ZERO
253C ENDDO
254C MODULE (EYN INITIALIZED TO E0)
255C DO J=10,12
256C UVAR(I,J)=E0
257C ENDDO
258C POISSON'S RATIO/MODULE
259C DO J=13,15
260C UVAR(I,J)=VT/E0
261C ENDDO
262C PRESSURE
263C UVAR(I,16)=ZERO
264C PRINCIPAL DIRECTIONS
265C DO J=17,25
266C UVAR(I,J)=ZERO
267C ENDDO
268C UVAR(I,17)=ONE
269C UVAR(I,21)=ONE
270C UVAR(I,25)=ONE
271C LAW STORAGE FOR DECAY AND HYSTERESIS
272C DO J=26,31
273C UVAR(I,J)=ZERO
274C ENDDO
275C ENDIF
276 ENDDO
277
278 IF(iflag==2)THEN
279 jacobi=.true.
280 iflag=0
281 ELSE
282 jacobi=.false.
283 ENDIF
284C......................................................................
285
286C DEFINE STRETCH TENSOR
287
288
289 IF(iflag==0)THEN
290C PRINCIPAL STRAINS FORMULATION
291 mdir=3
292 DO i=1,nel
293 IF(total) THEN
294 av(1,i)=epsxx(i)
295 av(2,i)=epsyy(i)
296 av(3,i)=epszz(i)
297 av(4,i)=epsxy(i)*half
298 av(5,i)=epsyz(i)*half
299 av(6,i)=epszx(i)*half
300 ELSE
301 av(1,i)=depsxx(i)
302 av(2,i)=depsyy(i)
303 av(3,i)=depszz(i)
304 av(4,i)=depsxy(i)*half
305 av(5,i)=depsyz(i)*half
306 av(6,i)=depszx(i)*half
307 ENDIF
308 ENDDO
309 IF(jacobi) THEN
310 DO i=1,nel
311 a(1,1)=av(1,i)
312 a(2,2)=av(2,i)
313 a(3,3)=av(3,i)
314 a(2,1)=av(4,i)
315 a(3,2)=av(5,i)
316 a(3,1)=av(6,i)
317 CALL jacobiew(a,3,earj,dpra,nrot)
318 earv(1,i)=earj(1)
319 earv(2,i)=earj(2)
320 earv(3,i)=earj(3)
321 DO n=1,3
322 DO j=1,3
323 dirprv(n,j,i)=dpra(n,j)
324 ENDDO
325 ENDDO
326 ENDDO
327 ELSE
328C Eigenvalues needed to be calculated in double precision
329C for a simple precision executing
330 IF (iresp==1) THEN
331 CALL valpvecdp(av,earv,dirprv,nel)
332 ELSE
333 CALL valpvec(av,earv,dirprv,nel)
334 ENDIF
335 ENDIF
336C SET OLD PRINCIPAL DIRECTIONS
337 DO i=1,nel
338 m=0
339 DO n=1,3
340 DO l=1,3
341 m=m+1
342 dprao(l,n)=uvar(i,16+m)
343 ENDDO
344 ENDDO
345 m=0
346 DO n=1,3
347 DO j=1,3
348 m=m+1
349 dpra(n,j)=dirprv(n,j,i)
350 ENDDO
351 ENDDO
352
353C COMPARE NEW AND OLD PRINCIPAL DIRECTIONS AND TRANSFORM STORED
354C VALUES FROM OLD INTO NEW PRINCIPAL DIRECTIONS IF TOLERANCE EXCEEDED
355
356 CALL checkaxes(dprao,dpra,amin,amax,unchanged,tolerance)
357 IF(.NOT.unchanged) THEN
358 DO m=1,4
359 DO n=1,3
360 DO l=1,3
361 sn(n,l)=zero
362 IF(n==l) sn(n,l)= uvar(i,n+3*(m-1))
363 ENDDO
364 ENDDO
365C ROTATE STORED VALUES IN OLD PRINCIPAL DIRECTIONS INTO GLOBAL DIRECTIONS
366 CALL dreh(sn,dprao,1,1,1)
367C ROTATE STORED VALUES IN GLOBAL DIRECTIONS INTO NEW PRINCIPAL DIRECTIONS
368 CALL dreh(sn,dpra ,1,1,0)
369 DO n=1,3
370 uvar(i,n+3*(m-1)) = sn(n,n)
371 ENDDO
372 ENDDO
373C pour la partie decay et hysteresis
374 IF (beta>tiny) THEN
375 DO m=1,2
376 DO n=1,3
377 DO l=1,3
378 sn(n,l)=zero
379 IF(n==l) sn(n,l)= uvar(i,n+3*(m-1)+25)
380 ENDDO
381 ENDDO
382 CALL dreh(sn,dprao ,1,1,1)
383 CALL dreh(sn,dpra ,1,1,0)
384 DO n=1,3
385 uvar(i,n+3*(m-1)+25) = sn(n,n)
386 ENDDO
387 ENDDO
388 ENDIF
389 uvar(i,17)=dpra(1,1)
390 uvar(i,18)=dpra(2,1)
391 uvar(i,19)=dpra(3,1)
392 uvar(i,20)=dpra(1,2)
393 uvar(i,21)=dpra(2,2)
394 uvar(i,22)=dpra(3,2)
395 uvar(i,23)=dpra(1,3)
396 uvar(i,24)=dpra(2,3)
397 uvar(i,25)=dpra(3,3)
398 ENDIF
399
400 IF(total) THEN
401 ear(i,1)=earv(1,i)
402 ear(i,2)=earv(2,i)
403 ear(i,3)=earv(3,i)
404 ELSE
405 ecr(i,1)=earv(1,i)
406 ecr(i,2)=earv(2,i)
407 ecr(i,3)=earv(3,i)
408 ENDIF
409 ENDDO
410 ELSE
411
412C AVERAGE STRAIN FORMULATION
413
414 mdir=1
415 DO i=1,nel
416 IF(total) THEN
417 ear(i,1)= sign(one,epsxx(i)+epsyy(i)+epszz(i))*sqrt(
418 & (epsxx(i)-epsyy(i))**two+
419 & (epsyy(i)-epszz(i))**two+
420 & (epszz(i)-epsxx(i))**two+
421 & two*(epsxy(i)**two+epsyz(i)**two+epszx(i)**two))/sqrt(two)
422 ELSE
423 ecr(i,1)= sign(one,depsxx(i)+depsyy(i)+depszz(i))*sqrt(
424 & (depsxx(i)-depsyy(i))**two+
425 & (depsyy(i)-depszz(i))**two+
426 & (depszz(i)-depsxx(i))**two+
427 & ((depsxy(i))**two+(depsyz(i))**two+(depszx(i))**two))
428 & /sqrt(two)
429 ENDIF
430 ENDDO
431
432 ENDIF
433
434C......................................................................
435 IF(timestep<=zero) THEN
436 dti=zero
437 ELSE
438 dti=one/timestep
439 ENDIF
440
441 DO j=1,mdir
442 DO i=1,nel
443 IF( total) ecr(i,j)=ear(i,j)-uvar(i,j)
444 IF(.NOT.total) ear(i,j)=ecr(i,j)+uvar(i,j)
445 ebr(i,j)=ecr(i,j)*dti
446 ENDDO
447C COMPUTE NOMINAL VALUES
448 DO i=1,nel
449C COMPUTE PRINCIPAL STRETCHES
450C LAMDA=E_NOMINAL+1=EXP(E_RADIOSS)
451 IF(ismstr==0.OR.ismstr==2.OR.ismstr==4) THEN
452 el(i,j)=exp(ear(i,j))
453C COMPUTE PRINCIPAL NOMINAL STRAIN(J) RATES FROM PRINCIPAL STRAIN(J) RATES
454C E_DOT_NOMINAL=EDOT_RADIOSS*EXP(E_RADIOSS)=EDOT_RADIOSS*LAMDA
455 ebn(i,j)=ebr(i,j)*el(i,j)
456 ecn(i,j)=el(i,j)*(one-exp(-ecr(i,j)))
457 ELSEIF(ismstr==10.OR.ismstr==12) THEN
458 el(i,j)=sqrt(ear(i,j)+one)
459 ebn(i,j)=ebr(i,j)
460 ecn(i,j)=ecr(i,j)
461 ELSE
462 el(i,j)=ear(i,j)+one
463 ebn(i,j)=ebr(i,j)
464 ecn(i,j)=ecr(i,j)
465 ENDIF
466C COMPUTE PRINCIPAL NOMINAL STRAINS FROM PRINCIPAL STRAINS
467 ean(i,j)=el(i,j)-one
468C STRAIN RATE FILTERING
469 ebn(i,j)=uvar(i,j+6)+(ebn(i,j)-uvar(i,j+6))*ratedamp
470 ENDDO
471 ENDDO
472
473C......................................................................
474
475C RELATIVE VOLUME
476 DO i=1,nel
477 pair(i)=zero
478 volumer(i)=el(i,1)*el(i,2)*el(i,3)
479 ENDDO
480
481C CONFINED AIR PRESSURE
482 DO i=1,nel
483 IF (kcompair==1.AND.volumer(i)<one) THEN
484 npcurve=ifunc(nfuncp)
485 IF (volumer(i)-phi>small) THEN
486 IF(npcurve/=0) THEN
487 pair(i)=-finter(npcurve,volumer(i),npf,tf,df)
488 pair(i)=pair(i)*pscale
489 ELSE
490 pair(i)=p0*(volumer(i)-one)/(volumer(i)-phi)
491 ENDIF
492 ELSE
493 pair(i)=uvar(i,16)
494 ENDIF
495C RELAXATION ON AIR PRESSURE
496 pair(i)=exp(-relaxp*time)*max(pair(i),-maxpres)
497 uvar(i,16)=pair(i)
498 ELSEIF (kcompair==2.AND.volumer(i)<one) THEN
499 npcurve=ifunc(nfuncp)
500 pair(i)=-finter(npcurve,volumer(i),npf,tf,df)
501 ENDIF
502 ENDDO
503
504C......................................................................
505
506C UPDATE TENSILE MODULE (INCREASES WITH COMPRESSION)
507
508C COMPUTE STRESS COMPONENTS VIA INTERPOLATION
509
510C LOOP ON PRINCIPAL DIRECTIONS
511 DO j=1,mdir
512 DO i=1,nel
513 psn(i,j)=zero
514C COMPRESSION IS POSITIVE!
515 strain(i,j)=-ean(i,j)
516 rate(i,j)=abs(ebn(i,j))
517 shift=zero
518 unloading=.false.
519 IF((ean(i,j)<zero.AND.ebn(i,j)>zero).OR.
520 & (ean(i,j)>zero.AND.ebn(i,j)<zero))unloading=.true.
521C COMPUTE ELASTIC MODULE FROM STATIC CURVE
522 psn1(i,j)=-alpha(1)*finter(ifunc(1),strain(i,j),npf,tf,df)
523C CURVE INTERPOLATION
524C EDOT0 = ACTUAL RATE
525C EDOTU = UNLOADING RATE
526C EDOTL = UPPER BOUND RATE
527C EDOTS = LOWERBOUND RATE
528 edot0(i,j)=rate(i,j)
529 IF(edot0(i,j)<edot(1)) edot0(i,j)=edot(1)
530 IF(unloading.AND.iunload/=0) THEN
531C CASE OF IUNLOAD/=0 (UNLOAD BETWEEN THE IUNLOAD CURVE (UPPER LIMIT)
532C AND THE FIRST CURVE (LOWER LIMIT)- STATIC- )
533 IF (abs(runload-edot(1))<em20) THEN
534 psn(i,j)=-funload*
535 & (finter(ifunc(nfuncul),strain(i,j),npf,tf,df1))
536 eyn(i,j)=funload*df1
537 visc(i,j)=zero
538 ELSE
539 edotu(i,j)=edot0(i,j)
540 edots(i,j)=edot(1)
541 edotl(i,j)=max(runload,edots(i,j))
542
543 IF(edotu(i,j)<edots(i,j))edotu(i,j)=edots(i,j)
544 IF(edotu(i,j)>edotl(i,j))edotu(i,j)=edotl(i,j)
545
546 psn1(i,j)=-alpha(1)*
547 & finter(ifunc(1),strain(i,j),npf,tf,df1)
548 eyn(i,j)=alpha(1)*df1
549 psn2(i,j)=-funload*
550 & (finter(ifunc(nfuncul),strain(i,j),npf,tf,df2))
551
552 psn(i,j)=psn2(i,j)+(psn1(i,j)-psn2(i,j))*
553 & (one-((min(one,(edotu(i,j)-edots(i,j))/
554 & (max(tiny,edotl(i,j)-edots(i,j)))))**exponas))**exponbs
555
556 IF(abs(edotu(i,j)-edots(i,j))>tiny.AND.
557 & abs(psn(i,j)-psn1(i,j))>tiny)visc(i,j)=
558 & abs((psn(i,j)-psn1(i,j))/(edotu(i,j)-edots(i,j)))
559
560 visc(i,j)=min(visc(i,j),viscosity)
561 psn(i,j)=psn1(i,j)-visc(i,j)*(edotu(i,j)-edots(i,j))
562 ENDIF
563
564 ELSE
565 kun(i,j)=0
566 IF(unloading) kun(i,j)=nfunc1
567C IF(UNLOADING) SHIFT=TWO*(SIGN(1,STRAIN(I,J))-STRAIN(I,J))
568C CASE OF ONLY ONE CURVE INPUT I.E. NO STRAIN(J) RATE EFFECT LOADING
569C AND UNLOADING FOLLOW THE SAME PATH
570 IF (nfunc1==1) THEN
571 psn(i,j)=-alpha(1)*
572 & finter(ifunc(kun(i,j)+1),strain(i,j)+shift,npf,tf,df)
573 eyn(i,j)=alpha(1)*df
574 ELSE
575C CASE OF IUNLOAD==0 ( UNLOAD BETWEEN THE IUNLOAD CURVE (UPPER LIMIT)
576C AND THE FIRST CURVE (LOWER LIMIT )- STATIC-
577C IF UNLOADING CURVES ARE SPECIFIED THEY WILL BE USED FOR INTERPOLATION
578C OTHERWISE STARTER INITIALIZES THEM TO THE DEFAULT CURVE NUMBER I.E.
579C THE FIRST LOADING (STATIC) CURVE
580 DO l=1,nfunc1
581 IF(edot0(i,j)<=edot(l)) GOTO 10
582 ENDDO
583 l=nfunc1
584 edot0(i,j)=edot(l)
585 10 CONTINUE
586 k(i,j)=l
587 kf(i,j)=l+kun(i,j)
588
589 psn2(i,j)=-alpha(k(i,j))*
590 & finter(ifunc(kf(i,j)),strain(i,j)+shift,npf,tf,df2)
591 k1(i,j)=max(1,l-1)
592 kf1(i,j)=max(1,l-1)+kun(i,j)
593 edotl(i,j)=edot(k(i,j))
594 edots(i,j)=edot(k1(i,j))
595 psn1(i,j)=-alpha(k1(i,j))*
596 & finter(ifunc(kf1(i,j)),strain(i,j)+shift,npf,tf,df1)
597 IF(l==nfunc1) THEN
598 eyn(i,j)=alpha(l)*df2
599 ELSE
600 eyn(i,j)=alpha(k1(i,j))*df1
601 ENDIF
602
603 psn(i,j)=psn2(i,j)+(psn1(i,j)-psn2(i,j))*
604 & (one-((min(one,(edot0(i,j)-edots(i,j))/
605 & (max(tiny,edotl(i,j)-edots(i,j)))))**exponas))
606 & **exponbs
607
608 IF(abs(edot0(i,j)-edots(i,j))>tiny.AND.
609 & abs(psn(i,j)-psn1(i,j))>tiny)visc(i,j)=
610 & abs((psn(i,j)-psn1(i,j))/(edot0(i,j)-edots(i,j)))
611
612 visc(i,j)=min(visc(i,j),viscosity)
613 psn(i,j)=psn1(i,j)-visc(i,j)*(edot0(i,j)-edots(i,j))
614 ENDIF
615 ENDIF
616
617 ENDDO
618
619C......................................................................
620
621C UNLOADING
622
623 DO i=1,nel
624 decay = uvar(i,j+28)
625C HYSTERESIS AND DECAY OCCURS ON COMPRESSION ONLY BUT REMAINS PERMANENT
626C IN TENSION
627 IF(ean(i,j)<zero) THEN
628C NO RECOVERY ON UNLOADING - ACCUMULATE STRAIN
629 IF(krecover==0.AND.ecn(i,j)<zero)
630 & uvar(i,j+25)=uvar(i,j+25)+abs(ecn(i,j))
631C RECOVERY ON UNLOADING - DOWN INCREMENT STRAIN
632 IF(krecover==1)uvar(i,j+25)=uvar(i,j+25)-ecn(i,j)
633 IF(ebn(i,j)<zero) THEN
634 decay = min(one,hyster*(one-exp(-beta*uvar(i,j+25))))
635 ENDIF
636 ENDIF
637C DECAY ON LOADING AND UNLOADING
638 IF(kdecay==0)
639 & psn(i,j)=psn(i,j)*(one-decay)
640C DECAY ON LOADING ONLY
641 IF(kdecay==1.AND.ecn(i,j)<zero)
642 & psn(i,j)=psn(i,j)*(one-decay)
643C DECAY ON UNLOADING ONLY
644 IF(kdecay==2.AND.ecn(i,j)>zero)
645 & psn(i,j)=psn(i,j)*(one-decay)
646 uvar(i,j+28)=decay
647C......................................................................
648 dsigma(i,j)=psn(i,j)-uvar(i,j+3)
649 IF(total) THEN
650C SECANT MODULE
651 tmp1 = eyn(i,j)
652 ELSE
653 IF(abs(ecn(i,j))>tiny)THEN
654 eyn(i,j)=abs(dsigma(i,j)/ecn(i,j))
655 ELSE
656 eyn(i,j)=uvar(i,j+9)
657 ENDIF
658 ENDIF
659 ENDDO
660
661 IF(.NOT.total) THEN
662 DO i=1,nel
663 IF(sign(one,(ebn(i,j)/(uvar(i,j+6)+tiny)))/=one)
664 & eyn(i,j)=uvar(i,j+9)
665 ENDDO
666
667 DO i=1,nel
668 eyn(i,j)=eyn(i,j)*theta+uvar(i,j+9)*(one-theta)
669 ENDDO
670 ENDIF
671
672
673 IF(itotal==0.OR.itotal==2) THEN
674C LINEAR IN TENSION
675 DO i=1,nel
676 IF(ean(i,j)>=zero) THEN
677 uvar(i,j+28)=zero
678 tmp1=exp(-lamda*(volumer(i)-1.+epsfin))
679 ei(i)=efinal+(e0-efinal)*(1-tmp1)
680 tmp2=lamda*(efinal-e0)*tmp1
681 eyn(i,j)=max(ei(i),tmp2)
682 IF(total) THEN
683 psn(i,j)=ei(i)*ean(i,j)
684 ELSE
685 psn(i,j)=uvar(i,j+3)+ei(i)*ecn(i,j)
686 ENDIF
687 uvar(i,j+12)=vt/ei(i)
688 ENDIF
689 ENDDO
690 ELSEIF(itotal==-1) THEN
691C MODULE MAXMAL IN TENSION ONE CYCLE LATER (TOTAL)
692 DO i=1,nel
693 IF(ean(i,j)>0.) THEN
694 uvar(i,j+28)=zero
695 ei(i)=max(uvar(i,10),uvar(i,11),uvar(i,12))
696 eyn(i,j)=ei(i)
697 psn(i,j)=ei(i)*ean(i,j)
698 uvar(i,j+12)=vt/ei(i)
699 ENDIF
700 ENDDO
701 ELSEIF(itotal==-2) THEN
702C MODULE MAXMAL IN TENSION ONE CYCLE LATER (INCR)
703 DO i=1,nel
704 IF(ean(i,j)>zero) THEN
705 uvar(i,j+28)=zero
706 ei(i)=max(uvar(i,10),uvar(i,11),uvar(i,12))
707 eyn(i,j)=ei(i)
708 psn(i,j)=uvar(i,j+3)+ei(i)*ecn(i,j)
709 uvar(i,j+12)=vt/ei(i)
710 ENDIF
711 ENDDO
712 ENDIF
713
714C END OF LOOP ON MDIR
715 ENDDO
716C BEGIN OF LOOP ON MDIR
717 IF (imsta>=1) THEN
718 DO i=1,nel
719 sigmax(i)=-(min(psn(i,1),psn(i,2),psn(i,3))-
720 . max(psn(i,1),psn(i,2),psn(i,3)))
721 ENDDO
722 DO j=1,mdir
723 DO i=1,nel
724 esecant(i,j)=0.4*abs(psn(i,j))/max(tiny,abs(strain(i,j)))
725 IF (esecant(i,j)<=sigmax(i)) THEN
726 tmp1=0.2*(sigmax(i)-esecant(i,j))
727 psn(i,j)=psn(i,j)+tmp1*ean(i,j)
728 eyn(i,j)=max(eyn(i,j),(one+tmp1)*esecant(i,j))
729 ENDIF
730 ENDDO
731 ENDDO
732 ENDIF
733 DO j=1,mdir
734 DO i=1,nel
735 tmp0=max(eyn(i,j),uvar(i,j+9))
736 uvar(i,j )=ear(i,j)
737 uvar(i,j+3)=psn(i,j)
738 uvar(i,j+6)=ebn(i,j)
739 uvar(i,j+9)=eyn(i,j)
740 eyn(i,j)=tmp0
741 ENDDO
742C END OF LOOP ON MDIR
743 ENDDO
744
745 IF(iflag==1) THEN
746C OCTAHEDRAL STRAIN (AVERAGE) FORMULATION
747 DO i=1,nel
748 emax(i)=ei(i)
749 d11(i)=emax(i)*(one-vt)/(one+vt)/(one-two*vt)
750 d12(i)=emax(i)*vt/(one+vt)/(one-two*vt)
751 d44(i)=emax(i)/two/(one+vt)
752 ENDDO
753 DO i=1,nel
754 IF(total) THEN
755 signxx(i)=d11(i)*epsxx(i)+d12(i)*epsyy(i)+d12(i)*epszz(i)
756 signyy(i)=d12(i)*epsxx(i)+d11(i)*epsyy(i)+d12(i)*epszz(i)
757 signzz(i)=d12(i)*epsxx(i)+d12(i)*epsyy(i)+d11(i)*epszz(i)
758 signxy(i)=d44(i)*epsxy(i)
759 signyz(i)=d44(i)*epsyz(i)
760 signzx(i)=d44(i)*epszx(i)
761 soundsp(i) = sqrt(d11(i)/rho0(i))
762 ELSE
763 signxx(i)=sigoxx(i)+
764 & d11(i)*depsxx(i)+d12(i)*depsyy(i)+d12(i)*depszz(i)
765 signyy(i)=sigoyy(i)+
766 & d12(i)*depsxx(i)+d11(i)*depsyy(i)+d12(i)*depszz(i)
767 signzz(i)=sigozz(i)+
768 & d12(i)*depsxx(i)+d12(i)*depsyy(i)+d11(i)*depszz(i)
769 signxy(i)=sigoxy(i)+d44(i)*depsxy(i)
770 signyz(i)=sigoyz(i)+d44(i)*depsyz(i)
771 signzx(i)=sigozx(i)+d44(i)*depszx(i)
772 soundsp(i) = sqrt(d11(i)/rho0(i))
773 ENDIF
774 viscmax(i)=visc(i,1)
775
776 ENDDO
777 RETURN
778 ENDIF
779
780 DO i=1,nel
781 IF(vt+vc<=two*tiny) THEN
782 IF(total) THEN
783 psc(i,1)=psn(i,1)
784 psc(i,2)=psn(i,2)
785 psc(i,3)=psn(i,3)
786 ELSE
787 psc(i,1)=eyn(i,1)*ecn(i,1)
788 psc(i,2)=eyn(i,2)*ecn(i,2)
789 psc(i,3)=eyn(i,3)*ecn(i,3)
790 ENDIF
791 ELSE
792 e12(i)=(ean(i,1)+ean(i,2))/two
793 e23(i)=(ean(i,2)+ean(i,3))/two
794 e31(i)=(ean(i,3)+ean(i,1))/two
795 v12(i)=vc+(vt-vc)*(one-exp(-rv*abs(e12(i))))*
796 & (sign(one,e12(i))+one)/two
797 v23(i)=vc+(vt-vc)*(one-exp(-rv*abs(e23(i))))*
798 & (sign(one,e23(i))+one)/two
799 v31(i)=vc+(vt-vc)*(one-exp(-rv*abs(e31(i))))*
800 & (sign(one,e31(i))+one)/two
801 IF(total) THEN
802 detc(i)=one-v23(i)*v23(i)-v31(i)*v31(i)-v12(i)*v12(i)-
803 & two*v12(i)*v31(i)*v23(i)
804 a11(i)=(one-v23(i)*v23(i))
805 a12(i)=(v12(i)+v23(i)*v31(i))
806 a13(i)=(v31(i)+v23(i)*v12(i))
807 a22(i)=(one-v31(i)*v31(i))
808 a23(i)=(v23(i)+v31(i)*v12(i))
809 a33(i)=(one-v12(i)*v12(i))
810 psc(i,1)=(a11(i)*psn(i,1)+a12(i)*psn(i,2)+a13(i)*psn(i,3))/
811 & detc(i)
812 psc(i,2)=(a12(i)*psn(i,1)+a22(i)*psn(i,2)+a23(i)*psn(i,3))/
813 & detc(i)
814 psc(i,3)=(a13(i)*psn(i,1)+a23(i)*psn(i,2)+a33(i)*psn(i,3))/
815 & detc(i)
816 ELSE
817 uvar(i,13)=theta*v23(i)/ei(i)+(one-theta)*uvar(i,13)
818 uvar(i,14)=theta*v31(i)/ei(i)+(one-theta)*uvar(i,14)
819 uvar(i,15)=theta*v12(i)/ei(i)+(one-theta)*uvar(i,15)
820 detc(i)=one/(eyn(i,1)*eyn(i,2)*eyn(i,3))-
821 & uvar(i,13)*uvar(i,13)/eyn(i,1)-
822 & uvar(i,14)*uvar(i,14)/eyn(i,2)-
823 & uvar(i,15)*uvar(i,15)/eyn(i,3)-
824 & 2*uvar(i,13)*uvar(i,14)*uvar(i,15)
825 a11(i)=one/(eyn(i,2)*eyn(i,3))-uvar(i,13)*uvar(i,13)
826 a12(i)=uvar(i,15)/eyn(i,3)+uvar(i,13)*uvar(i,14)
827 a13(i)=uvar(i,14)/eyn(i,2)+uvar(i,13)*uvar(i,15)
828 a22(i)=one/(eyn(i,1)*eyn(i,3))-uvar(i,14)*uvar(i,14)
829 a23(i)=uvar(i,13)/eyn(i,1)+uvar(i,14)*uvar(i,15)
830 a33(i)=one/(eyn(i,1)*eyn(i,2))-uvar(i,15)*uvar(i,15)
831C COMPUTE STRESS INCREMENT IN PRINCIPAL DIRECTIONS
832 psc(i,1)=((a11(i)*ecn(i,1)+a12(i)*ecn(i,2)+a13(i)*ecn(i,3))
833 & /detc(i))
834 psc(i,2)=((a12(i)*ecn(i,1)+a22(i)*ecn(i,2)+a23(i)*ecn(i,3))
835 & /detc(i))
836 psc(i,3)=((a13(i)*ecn(i,1)+a23(i)*ecn(i,2)+a33(i)*ecn(i,3))
837 & /detc(i))
838 ENDIF
839 ENDIF
840
841 DO j=1,3
842 IF(off(i)==zero.OR.psn(i,j)>abs(tensioncut))THEN
843 psc(i,1)=zero
844 psc(i,2)=zero
845 psc(i,3)=zero
846 off(i)=zero
847 ENDIF
848 ENDDO
849C COMPUTE CAUCHY STRESS INCREMENT IN PRINCIPAL DIRECTIONS
850CSIGMA_CAUCHY(I) = LAMDA(I) * SIGMA_NOMINAL(I) / RELATIVE VOLUME
851 IF(ismstr==0.OR.ismstr==2.OR.ismstr==4) THEN
852C CAUCHY MODULE FOR TIME STEP
853 den1=el(i,2)*el(i,3)
854 den2=el(i,3)*el(i,1)
855 den3=el(i,1)*el(i,2)
856 psc(i,1)=psc(i,1)/den1
857 psc(i,2)=psc(i,2)/den2
858 psc(i,3)=psc(i,3)/den3
859 eyn(i,1)=eyn(i,1)/den1
860 eyn(i,2)=eyn(i,2)/den2
861 eyn(i,3)=eyn(i,3)/den3
862 ELSEIF(ismstr==10.OR.ismstr==12) THEN
863C------------directly to Chauchy-------
864 den1=el(i,1)/volumer(i)
865 den2=el(i,2)/volumer(i)
866 den3=el(i,3)/volumer(i)
867 psc(i,1)=psc(i,1)*den1
868 psc(i,2)=psc(i,2)*den2
869 psc(i,3)=psc(i,3)*den3
870 eyn(i,1)=eyn(i,1)*den1
871 eyn(i,2)=eyn(i,2)*den2
872 eyn(i,3)=eyn(i,3)*den3
873 ENDIF
874 ENDDO
875 IF (kcompair==2) THEN
876 DO i=1,nel
877 tmp0=volumer(i)
878 tmp3=min(el(i,1),el(i,2),el(i,3))
879 IF (tmp0<one.AND.tmp3<one
880 & .AND.tmp3>tmp0.AND.abs(tmp0-tmp3)>em6) THEN
881 tmp2=volumer(i)**(1./3.)-volumer(i)
882 tmp4=(tmp3-tmp0)/tmp2
883 aa(i)=min(one,tmp4)
884 ELSE
885 aa(i)=zero
886 ENDIF
887 ENDDO
888 DO j=1,3
889 DO i=1,nel
890 sigprv(j,i)=psc(i,j)+aa(i)*(pair(i)-psc(i,j))
891 ENDDO
892 ENDDO
893 ELSE
894 DO j=1,3
895 DO i=1,nel
896 sigprv(j,i)=psc(i,j)+pair(i)
897 ENDDO
898 ENDDO
899 ENDIF
900C TRANSFORM FROM PRINCIPAL TO GLOBAL DIRECTIONS
901 DO i=1,nel
902 signxx(i) = dirprv(1,1,i)*dirprv(1,1,i)*sigprv(1,i)
903 & + dirprv(1,2,i)*dirprv(1,2,i)*sigprv(2,i)
904 & + dirprv(1,3,i)*dirprv(1,3,i)*sigprv(3,i)
905 signyy(i) = dirprv(2,2,i)*dirprv(2,2,i)*sigprv(2,i)
906 & + dirprv(2,3,i)*dirprv(2,3,i)*sigprv(3,i)
907 & + dirprv(2,1,i)*dirprv(2,1,i)*sigprv(1,i)
908 signzz(i) = dirprv(3,3,i)*dirprv(3,3,i)*sigprv(3,i)
909 & + dirprv(3,1,i)*dirprv(3,1,i)*sigprv(1,i)
910 & + dirprv(3,2,i)*dirprv(3,2,i)*sigprv(2,i)
911 signxy(i) = dirprv(1,1,i)*dirprv(2,1,i)*sigprv(1,i)
912 & + dirprv(1,2,i)*dirprv(2,2,i)*sigprv(2,i)
913 & + dirprv(1,3,i)*dirprv(2,3,i)*sigprv(3,i)
914 signyz(i) = dirprv(2,2,i)*dirprv(3,2,i)*sigprv(2,i)
915 & + dirprv(2,3,i)*dirprv(3,3,i)*sigprv(3,i)
916 & + dirprv(2,1,i)*dirprv(3,1,i)*sigprv(1,i)
917 signzx(i) = dirprv(3,3,i)*dirprv(1,3,i)*sigprv(3,i)
918 & + dirprv(3,1,i)*dirprv(1,1,i)*sigprv(1,i)
919 & + dirprv(3,2,i)*dirprv(1,2,i)*sigprv(2,i)
920 ENDDO
921
922 IF(.NOT.total) THEN
923 DO i=1,nel
924C ADD INCREMENT OF STRESS
925 signxx(i)=signxx(i)+sigoxx(i)
926 signyy(i)=signyy(i)+sigoyy(i)
927 signzz(i)=signzz(i)+sigozz(i)
928 signxy(i)=signxy(i)+sigoxy(i)
929 signyz(i)=signyz(i)+sigoyz(i)
930 signzx(i)=signzx(i)+sigozx(i)
931 ENDDO
932 ENDIF
933C SOUNDSPEED
934 IF(iflag==0) THEN
935 DO i=1,nel
936 emax(i)=max(ei(i),eyn(i,1),eyn(i,2),eyn(i,3))
937 ENDDO
938 ENDIF
939
940 IF (imsta==2) THEN
941 DO i=1,nel
942 epsxy(i) = dirprv(1,1,i)*dirprv(2,1,i)*ean(i,1)
943 & + dirprv(1,2,i)*dirprv(2,2,i)*ean(i,2)
944 & + dirprv(1,3,i)*dirprv(2,3,i)*ean(i,3)
945 epsyz(i) = dirprv(2,2,i)*dirprv(3,2,i)*ean(i,2)
946 & + dirprv(2,3,i)*dirprv(3,3,i)*ean(i,3)
947 & + dirprv(2,1,i)*dirprv(3,1,i)*ean(i,1)
948 epszx(i) = dirprv(3,3,i)*dirprv(1,3,i)*ean(i,3)
949 & + dirprv(3,1,i)*dirprv(1,1,i)*ean(i,1)
950 & + dirprv(3,2,i)*dirprv(1,2,i)*ean(i,2)
951 ENDDO
952 DO i=1,nel
953 esecant(i,1)=0.5*abs(signxy(i))/max(tiny,abs(epsxy(i)))
954 esecant(i,2)=0.5*abs(signyz(i))/max(tiny,abs(epsyz(i)))
955 esecant(i,3)=0.5*abs(signzx(i))/max(tiny,abs(epszx(i)))
956 sigmax(i)=max(0.5*ei(i),sigmax(i))
957 IF (esecant(i,1)<=sigmax(i)) THEN
958 tmp1=0.1*(sigmax(i)-esecant(i,1))
959 signxy(i)=signxy(i)+tmp1*epsxy(i)
960 ENDIF
961 IF (esecant(i,2)<=sigmax(i)) THEN
962 tmp1=0.1*(sigmax(i)-esecant(i,2))
963 signyz(i)=signyz(i)+tmp1*epsyz(i)
964 ENDIF
965 IF (esecant(i,3)<=sigmax(i)) THEN
966 tmp1=0.1*(sigmax(i)-esecant(i,3))
967 signzx(i)=signzx(i)+tmp1*epszx(i)
968 ENDIF
969 ENDDO
970 ENDIF
971 DO i=1,nel
972 kkk=emax(i)/(three*(one-two*max(vc,vt)))
973 ggg=emax(i)/(two*(one+max(vc,vt)))
974 soundsp(i)=sqrt((kkk+four_over_3*ggg)/rho0(i))
975 viscmax(i)=max(visc(i,1),visc(i,2),visc(i,3))
976 ENDDO
977
978 RETURN
979
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
valpvecdp
subroutine valpvecdp(sig, val, vec, nel)
Definition matrix.F:533
jacobiew
subroutine jacobiew(a, n, ew, ev, nrot)
Definition matrix.F:29
valpvec
subroutine valpvec(sig, val, vec, nel)
Definition matrix.F:215
checkaxes
subroutine checkaxes(tran1, tran2, amin, amax, unchanged, tol)
Definition sigeps38.F:988