OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps90.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!|| sigeps90 ../engine/source/materials/mat/mat090/sigeps90.F
25!||--- called by ------------------------------------------------------
26!|| mulaw ../engine/source/materials/mat_share/mulaw.F90
27!||--- calls -----------------------------------------------------
28!|| valpvec_v ../engine/source/materials/mat/mat033/sigeps33.F
29!|| valpvecdp_v ../engine/source/materials/mat/mat033/sigeps33.F
30!|| vinter2 ../engine/source/tools/curve/vinter.F
31!||--- uses -----------------------------------------------------
32!|| file_descriptor_mod ../engine/source/modules/file_descriptor_mod.F90
33!||====================================================================
34 SUBROUTINE sigeps90(
35 1 NEL ,NUVAR ,NFUNC ,IFUNC ,NPF ,
36 2 TF ,TIME ,UPARAM ,RHO0 ,
37 3 EPSPXX ,EPSPYY ,EPSPZZ ,EPSPXY ,EPSPYZ ,EPSPZX ,
38 4 EPSXX ,EPSYY ,EPSZZ ,EPSXY ,EPSYZ ,EPSZX ,
39 5 SIGOXX ,SIGOYY ,SIGOZZ ,SIGOXY ,SIGOYZ ,SIGOZX ,
40 5 SIGNXX ,SIGNYY ,SIGNZZ ,SIGNXY ,SIGNYZ ,SIGNZX ,
41 6 SOUNDSP,VISCMAX,UVAR ,NVARTMP ,VARTMP ,ISMSTR ,
42 7 ISRATE ,ASRATE ,OFFG ,IHET ,ET ,EPSD ,
43 8 OFF , NGL )
44! ----------------------------------------------------------------------------------------------------------------------
45! Modules
46! ----------------------------------------------------------------------------------------------------------------------
47 use file_descriptor_mod
48C-----------------------------------------------
49C I m p l i c i t T y p e s
50C-----------------------------------------------
51#include "implicit_f.inc"
52#include "comlock.inc"
53#include "mvsiz_p.inc"
54C-----------------------------------------------
55C G l o b a l P a r a m e t e r s
56C-----------------------------------------------
57C---------+---------+---+---+--------------------------------------------
58C VAR | SIZE |TYP| RW| DEFINITION
59C---------+---------+---+---+--------------------------------------------
60C NEL | 1 | I | R | SIZE OF THE ELEMENT GROUP NEL
61C NUPARAM | 1 | I | R | SIZE OF THE USER PARAMETER ARRAY
62C NUVAR | 1 | I | R | NUMBER OF USER ELEMENT VARIABLES
63C---------+---------+---+---+--------------------------------------------
64C MFUNC | 1 | I | R | NUMBER FUNCTION USED FOR THIS USER LAW not used
65C KFUNC | NFUNC | I | R | FUNCTION INDEX not used
66C NPF | * | I | R | FUNCTION ARRAY
67C TF | * | F | R | FUNCTION ARRAY
68C---------+---------+---+---+--------------------------------------------
69C TIME | 1 | F | R | CURRENT TIME
70C UPARAM | NUPARAM | F | R | USER MATERIAL PARAMETER ARRAY
71C RHO0 | NEL | F | R | INITIAL DENSITY
72C RHO | NEL | F | R | DENSITY
73C EINT | NEL | F | R | TOTAL INTERNAL ENERGY
74C EPSPXX | NEL | F | R | STRAIN RATE XX
75C EPSPYY | NEL | F | R | STRAIN RATE YY
76C ... | | | |
77C EPSXX | NEL | F | R | STRAIN XX
78C EPSYY | NEL | F | R | STRAIN YY
79C ... | | | |
80C---------+---------+---+---+--------------------------------------------
81C SIGNXX | NEL | F | W | NEW ELASTO PLASTIC STRESS XX
82C SIGNYY | NEL | F | W | NEW ELASTO PLASTIC STRESS YY
83C ... | | | |
84C SIGVXX | NEL | F | W | VISCOUS STRESS XX
85C SIGVYY | NEL | F | W | VISCOUS STRESS YY
86C ... | | | |
87C SOUNDSP | NEL | F | W | SOUND SPEED (NEEDED FOR TIME STEP)
88C VISCMAX | NEL | F | W | MAXIMUM DAMPING MODULUS(NEEDED FOR TIME STEP)
89C---------+---------+---+---+--------------------------------------------
90C UVAR |NEL*NUVAR| F |R/W| USER ELEMENT VARIABLE ARRAY
91C OFF | NEL | F |R/W| DELETED ELEMENT FLAG (=1. ON, =0. OFF)
92C---------+---------+---+---+--------------------------------------------
93#include "scr05_c.inc"
94#include "impl1_c.inc"
95C
96 INTEGER, INTENT(IN) :: NEL,NUVAR,ISMSTR,ISRATE,IHET
97 INTEGER, INTENT(IN) :: NVARTMP,NGL(NEL)
98
99 my_real, INTENT(IN) ,DIMENSION(NEL) :: RHO0,EPSPXX,EPSPYY,EPSPZZ,
100 . EPSPXY,EPSPYZ,EPSPZX,
101 . EPSXX,EPSYY,EPSZZ,
102 . EPSXY ,EPSYZ,EPSZX ,
103 . OFFG,SIGOXX, SIGOYY,
104 . SIGOZZ, SIGOXY, SIGOYZ,SIGOZX
105 my_real, INTENT(IN) :: time,uparam(*),asrate
106C-----------------------------------------------
107C O U T P U T A r g u m e n t s
108C-----------------------------------------------
109 my_real,INTENT(OUT)::
110 . signxx(nel),signyy(nel),signzz(nel),
111 . signxy(nel),signyz(nel),signzx(nel),
112 . soundsp(nel),viscmax(nel),et(nel),epsd(nel),off(nel)
113C-----------------------------------------------
114C I N P U T O U T P U T A r g u m e n t s
115C-----------------------------------------------
116 my_real, INTENT(INOUT) :: uvar(nel,nuvar)
117 INTEGER, DIMENSION(NEL,NVARTMP), INTENT(INOUT) :: VARTMP ! last interpolation positions in function tables
118C-----------------------------------------------
119C VARIABLES FOR FUNCTION INTERPOLATION
120C-----------------------------------------------
121 INTEGER, INTENT(IN) :: NPF(*), NFUNC, IFUNC(NFUNC)
122 my_real, INTENT(IN) :: TF(*)
123C-----------------------------------------------
124C L o c a l V a r i a b l e s
125C-----------------------------------------------
126 INTEGER :: I,J,K,II,KK,I1,J1,J2,IFLAG,IDAM,NE_L,NE_UNL,TFLAG,NINDX,
127 . NINDX_PRINT,N,FAIL
128 INTEGER :: ILOAD(MVSIZ,3),INDX_L(MVSIZ),INDX_UNL(MVSIZ),INDX(MVSIZ),
129 . INDX_PRINT(NEL)
130 INTEGER ,DIMENSION(NEL) :: IAD,ILEN,IPOS,JJ,UNLOAD
131 my_real
132 . e0,aa,g,nu,shape,hys,
133 . yfac1,yfacj1,yfacj2,ratej1,ratej2, epse,ep1,
134 . ep2,ep3,ep4,ep5,ep6,ert11,ert12,ert13,ert21,
135 . ert22,ert23,ert31,ert32,ert33,sj1,sj2,fac,t1,t2,t3,
136 . dam,epe,e_min(mvsiz),delta,alpha,tcut,de,rateeps,
137 . deps,e_max,e_new,e_old,epss,tcut0
138 my_real :: df(3),epsp(3)
139
140 my_real ,DIMENSION(NEL) :: dydx,stmp1,stmp2
141 my_real, DIMENSION(NEL,3) :: strain,strainrate,s,sqstat
142 my_real, DIMENSION(NEL,2) :: rate,yfac
143 my_real, DIMENSION(MVSIZ) :: quasi_eint,emax,emin,
144 . ecurent,e,deint,
145 . yld(mvsiz),epst(mvsiz)
146 my_real, DIMENSION(MVSIZ,6) :: av
147 my_real, DIMENSION(MVSIZ,6) :: evv,ev
148 my_real, DIMENSION(MVSIZ,3,3) :: dirprv
149C-----------------------------------------------
150C USER VARIABLES INITIALIZATION
151C-----------------------------------------------
152 e0 = uparam(1)
153 g = uparam(4)
154 nu = uparam(5)
155 shape = uparam(6)
156 hys = uparam(7)
157 iflag = uparam(9)
158 idam = uparam(10)
159 e_max = uparam(2*nfunc + 11)
160 alpha = uparam(2*nfunc + 13)
161 tflag = uparam(2*nfunc + 14)
162 fail = nint(uparam(2*nfunc + 15))
163 tcut0 = uparam(2*nfunc + 16)
164C
165 IF(time == zero )THEN
166 uvar(1:nel,8) = e0
167 uvar(1:nel,7) = one
168 ENDIF
169C
170C-----------------------------------------------
171C
172 ! Deletion of element
173 nindx = 0
174 indx(1:nel) = 0
175 nindx_print = 0
176 indx_print(1:nel) = 0
177 DO i=1,nel
178 IF(off(i) == zero ) THEN
179 signxx(i) = zero
180 signyy(i) = zero
181 signzz(i) = zero
182 signxy(i) = zero
183 signyz(i) = zero
184 signzx(i) = zero
185 e_old = uvar(i,8)
186 aa = e_old*(one-nu)/(one + nu)/(one - two*nu)
187 soundsp(i) = sqrt(aa/rho0(i))
188 ELSEIF(off(i) < one ) THEN
189 off(i) = off(i)*four_over_5
190 signxx(i) = sigoxx(i)*off(i)
191 signyy(i) = sigoyy(i)*off(i)
192 signzz(i) = sigozz(i)*off(i)
193 signxy(i) = sigoxy(i)*off(i)
194 signyz(i) = sigoyz(i)*off(i)
195 signzx(i) = sigozx(i)*off(i)
196 IF(off(i) < em01) THEN
197 off(i) = zero
198 signxx(i) = zero
199 signyy(i) = zero
200 signzz(i) = zero
201 signxy(i) = zero
202 signyz(i) = zero
203 signzx(i) = zero
204 e_old = uvar(i,8)
205 aa = e_old*(one-nu)/(one + nu)/(one - two*nu)
206 soundsp(i) = sqrt(aa/rho0(i))
207 ! for print_out
208 nindx_print = nindx_print + 1
209 indx_print(nindx_print) = i
210 ENDIF
211 ELSE
212 nindx = nindx + 1
213 indx(nindx) = i
214 ENDIF
215 ENDDO
216 !
217 IF (nindx_print > 0) THEN
218 DO j=1,nindx_print
219#include "lockon.inc"
220 WRITE(iout, 1000) ngl(indx_print(j))
221 WRITE(istdo,1100) ngl(indx_print(j)),time
222#include "lockoff.inc"
223 ENDDO
224 ENDIF
225 IF(nindx == 0 ) RETURN ! all groupe elemenet is deleted
226 !
227 ! Treatment of not deleted element.
228#include "vectorize.inc"
229 DO n=1,nindx
230 i = indx(n)
231 av(n,1) = epsxx(i)
232 av(n,2) = epsyy(i)
233 av(n,3) = epszz(i)
234 av(n,4) = half*epsxy(i)
235 av(n,5) = half*epsyz(i)
236 av(n,6) = half*epszx(i)
237 ENDDO
238C Eigenvalues needed to be calculated in double precision
239C for a simple precision executing*
240 IF (iresp==1) THEN
241 CALL valpvecdp_v(av,evv,dirprv,nindx)
242 ELSE
243 CALL valpvec_v(av,evv,dirprv,nindx)
244 ENDIF
245C-ISMSTR=0-NO SMALL STRAIN OPTION:STRAINS ARE LOGARITHMIC, STRESS IS CAUCHY
246C-ISMSTR=1-SMALL STRAIN OPTION:STRAINS ARE ENGINEERING, STRESS IS CAUCHY
247C-ISMSTR=2-SMALL STRAIN OPTION:STRAINS ARE ENGINEERING, STRESS IS BIOT
248C-ISMSTR=3-NO SMALL STRAIN OPTION:STRESS IS BIOT
249 IF(ismstr==0.OR.ismstr==2.OR.ismstr==4) THEN
250#include "vectorize.inc"
251 DO n=1,nindx
252 i = indx(n)
253C ---- (STRAIN IS LOGARITHMIC)
254 ev(n,1)=exp(evv(n,1))
255 ev(n,2)=exp(evv(n,2))
256 ev(n,3)=exp(evv(n,3))
257 ENDDO
258 ELSEIF(ismstr==10.OR.ismstr==12) THEN
259 DO n=1,nindx
260 i = indx(n)
261 IF(offg(i)<=one) THEN
262 ev(n,1)=sqrt(evv(n,1) + one )
263 ev(n,2)=sqrt(evv(n,2) + one )
264 ev(n,3)=sqrt(evv(n,3) + one )
265 ELSE
266 ev(n,1)=evv(n,1)+ one
267 ev(n,2)=evv(n,2)+ one
268 ev(n,3)=evv(n,3)+ one
269 END IF
270 ENDDO
271 ELSE
272C ---- STRAIN IS ENGINEERING)
273#include "vectorize.inc"
274 DO n=1,nindx
275 i = indx(n)
276 ev(n,1)=evv(n,1) + one
277 ev(n,2)=evv(n,2) + one
278 ev(n,3)=evv(n,3) + one
279 ENDDO
280 ENDIF
281C engineering strain and strain rate
282#include "vectorize.inc"
283 DO i=1,nindx
284 ii = indx(i)
285C engineering strain e = lambda-1 , according the input curve
286C e=1-lambda (e > 0 compression and e < 0 traction)
287 strain(i,1) = one - ev(i,1)
288 strain(i,2) = one - ev(i,2)
289 strain(i,3) = one - ev(i,3)
290
291 epst(i) = sqrt(strain(i,1)**2 + strain(i,2)**2 + strain(i,3)**2)
292C
293 ep1 = epspxx(ii)
294 ep2 = epspyy(ii)
295 ep3 = epspzz(ii)
296 ep4 = half*epspxy(ii)
297 ep5 = half*epspyz(ii)
298 ep6 = half*epspzx(ii)
299C phi_trans*L*phi_t
300 ert11 =dirprv(i,1,1)*ep1 + dirprv(i,2,1)*ep4 + dirprv(i,3,1)*ep6
301 ert12 =dirprv(i,1,2)*ep1 + dirprv(i,2,2)*ep4 + dirprv(i,3,2)*ep6
302 ert13 =dirprv(i,1,3)*ep1 + dirprv(i,2,3)*ep4 + dirprv(i,3,3)*ep6
303
304 ert21 =dirprv(i,1,1)*ep4 + dirprv(i,2,1)*ep2 + dirprv(i,3,1)*ep5
305 ert22 =dirprv(i,1,2)*ep4 + dirprv(i,2,2)*ep2 + dirprv(i,3,2)*ep5
306 ert23 =dirprv(i,1,3)*ep4 + dirprv(i,2,3)*ep2 + dirprv(i,3,3)*ep5
307
308 ert31 =dirprv(i,1,1)*ep6 + dirprv(i,2,1)*ep5 + dirprv(i,3,1)*ep3
309 ert32 =dirprv(i,1,2)*ep6 + dirprv(i,2,2)*ep5 + dirprv(i,3,2)*ep3
310 ert33 =dirprv(i,1,3)*ep6 + dirprv(i,2,3)*ep5 + dirprv(i,3,3)*ep3
311C
312 epsp(1) = dirprv(i,1,1)*ert11 + dirprv(i,2,1)*ert21
313 . + dirprv(i,3,1)*ert31
314 epsp(2) = dirprv(i,1,2)*ert12 + dirprv(i,2,2)*ert22
315 . + dirprv(i,3,2)*ert32
316 epsp(3) = dirprv(i,1,3)*ert13 + dirprv(i,2,3)*ert23
317 . + dirprv(i,3,3)*ert33
318C abs(eps) not necessary
319 strainrate(i,1) = epsp(1)*(one - strain(i,1)) ! eng
320 strainrate(i,2) = epsp(2)*(one - strain(i,2))
321 strainrate(i,3) = epsp(3)*(one - strain(i,3))
322 ENDDO
323!
324 ! computing energy increase
325 yfac1 = uparam(nfunc + 11)
326 quasi_eint(1:nindx)= zero
327!
328 DO k=1,3
329 kk = (k-1)*nfunc + 1
330 DO i=1,nindx
331 ii = indx(i)
332 ipos(i) = vartmp(ii,kk)
333 iad(i) = npf(ifunc(1)) / 2 + 1
334 ilen(i) = npf(ifunc(1)+1) / 2 - iad(i) - ipos(i)
335 ENDDO
336 CALL vinter2(tf,iad,ipos,ilen,nindx,strain(1,k),dydx,sqstat(1,k))
337 !
338 DO i=1,nindx
339 ii = indx(i)
340 vartmp(ii,kk) = ipos(i)
341 sqstat(i,k) = yfac1 * sqstat(i,k)
342 ENDDO
343 ! compute current energy
344 DO i=1,nindx
345 IF (tflag == 2 .AND. strain(i,k) < zero ) sqstat(i,k) = e0*strain(i,k)
346 quasi_eint(i) = quasi_eint(i) + half*strain(i,k)*sqstat(i,k)
347 ENDDO
348 ENDDO
349 !
350 DO i=1,nindx
351 ii = indx(i)
352 deint(i) = quasi_eint(i) - uvar(ii,9) !
353 uvar(ii,9) = quasi_eint(i)
354 ENDDO
355C -----------
356C check loading and unloading.
357 IF(nindx > 0 ) THEN
358 DO k=1,3
359 DO i=1,nindx
360 epe = epsp(k)*strain(i,k)
361 iload(i,k) = 1
362 IF(epe > em10 )iload(i,k) = -1
363 ENDDO
364 ENDDO
365 ENDIF
366C filtering strain rate
367 DO i=1,nindx
368 ii = indx(i)
369 rateeps = sqrt(strainrate(i,1)**2 + strainrate(i,2)**2 + strainrate(i,3)**2)
370 IF (israte > 0) THEN
371 rateeps = asrate*rateeps + (one - asrate)*uvar(ii,3)
372 ENDIF
373 epsd(ii) = rateeps
374 ENDDO
375C sous groupe
376 indx_l(1:nindx) = 0
377 indx_unl(1:nindx) = 0
378 ne_l = 0
379 ne_unl = 0
380 unload(1:nel) = 0
381C
382 DO i=1,nindx
383 ii = indx(i)
384 deps = epst(i) - uvar(ii,6)
385
386 IF(deint(i) >= zero .OR. deps >= zero) THEN
387 ne_l = ne_l + 1
388 indx_l(ne_l) = i
389 uvar(ii,3) = epsd(ii)
390 ELSE
391 epsd(ii) = min(epsd(ii), uvar(ii,3))
392 ne_unl = ne_unl + 1
393 indx_unl(ne_unl) = i
394 uvar(ii,3) = epsd(ii)
395 unload(i) = 1
396 ENDIF
397 strainrate(i,1) = epsd(ii)
398 strainrate(i,2) = epsd(ii)
399 strainrate(i,3) = epsd(ii)
400 ENDDO
401c
402C case with unloading stress-strain curve
403c
404 IF (iflag == 1) THEN
405 IF (nfunc == 1) THEN
406 yfac1 = uparam(nfunc + 11)
407 DO i=1,nel
408 emax(i) = zero
409 et(i) = zero
410 emin(i) = ep20
411 ENDDO
412!
413 DO k=1,3
414 kk = (k-1)*nfunc + 1
415 DO i=1,nindx
416 ii = indx(i)
417 ipos(i) = vartmp(ii,kk)
418 iad(i) = npf(ifunc(1)) / 2 + 1
419 ilen(i) = npf(ifunc(1)+1) / 2 - iad(i) - ipos(i)
420 ENDDO
421 CALL vinter2(tf,iad,ipos,ilen,nindx,strain(1,k),dydx,s(1,k))
422 DO i=1,nindx
423 ii = indx(i)
424 s(i,k) = yfac1 * s(i,k)
425 IF (tflag == 2 .AND. strain(i,k) < zero ) s(i,k) = e0*strain(i,k)
426 emax(i) = max(emax(i), yfac1*dydx(i))
427 emin(i) = max(emin(i), yfac1*dydx(i))
428 vartmp(ii,kk) = ipos(i)
429 ENDDO
430 ENDDO
431 DO i=1,nindx
432 ii = indx(i)
433 e(i) = emax(i)
434 et(ii) = emin(i)/e0
435 yld(i) = sqrt(s(i,1)**2 + s(i,2)**2 + s(i,3)**2)
436 ENDDO
437 ELSE ! multiple functions (strain rate interpolation
438 DO i=1,nindx
439 emax(i) = zero
440 emin(i) = ep20
441 ENDDO
442!
443 DO k=1,3 ! by direction
444 jj(1:nindx) = 1
445 DO i=1,nindx
446 DO j = 2,nfunc - 1
447 IF (strainrate(i,k) >= uparam(10 + j)) jj(i) = j
448 ENDDO
449 END DO
450 DO i=1,nindx
451 rate(i,1) = uparam(10 + jj(i))
452 rate(i,2) = uparam(10 + jj(i)+1)
453 yfac(i,1) = uparam(10 + nfunc + jj(i))
454 yfac(i,2) = uparam(10 + nfunc + jj(i)+1)
455 END DO
456 !
457 DO i=1,nindx
458 ii = indx(i)
459 j1 = jj(i)
460 kk = (k-1)*nfunc + j1
461 ipos(i) = vartmp(ii,kk)
462 iad(i) = npf(ifunc(j1)) / 2 + 1
463 ilen(i) = npf(ifunc(j1)+1) / 2 - iad(i) - ipos(i)
464 END DO
465 CALL vinter2(tf,iad,ipos,ilen,nindx,strain(1,k),dydx,stmp1)
466
467 DO i=1,nindx
468 ii = indx(i)
469 kk = (k-1)*nfunc + jj(i)
470 stmp1(i) = yfac(i,1) * stmp1(i)
471 emax(i) = max(emax(i), yfac(i,1)*dydx(i))
472 emin(i) = max(emin(i), yfac(i,1)*dydx(i))
473 vartmp(ii,kk) = ipos(i)
474 END DO
475
476 DO i=1,nindx
477 ii = indx(i)
478 j2 = jj(i)+1
479 kk = (k-1)*nfunc + j2
480 ipos(i) = vartmp(ii,kk)
481 iad(i) = npf(ifunc(j2)) / 2 + 1
482 ilen(i) = npf(ifunc(j2)+1) / 2 - iad(i) - ipos(i)
483 END DO
484 CALL vinter2(tf,iad,ipos,ilen,nindx,strain(1,k),dydx,stmp2)
485
486 DO i=1,nindx
487 ii = indx(i)
488 kk = (k-1)*nfunc + jj(i)+1
489 stmp2(i) = yfac(i,2) * stmp2(i)
490 emax(i) = max(emax(i), yfac(i,2)*dydx(i))
491 emin(i) = max(emin(i), yfac(i,2)*dydx(i))
492 vartmp(ii,kk) = ipos(i)
493 END DO
494!
495 DO i=1,nindx
496 fac = (strainrate(i,k) - rate(i,1)) / (rate(i,2) - rate(i,1))
497 s(i,k) = stmp1(i) + fac*(stmp2(i) - stmp1(i))
498 IF (tflag == 2 .AND. strain(i,k) < zero ) s(i,k) = e0*strain(i,k)
499 END DO
500
501 END DO ! K=1,3
502 DO i=1,nindx
503 ii = indx(i)
504 e(i) = max(e0,emax(i))
505 e(i) = emax(i)
506 et(ii) = emin(i)/e0 ! Not used
507 yld(i) = sqrt(s(i,1)**2 + s(i,2)**2 + s(i,3)**2)
508 ENDDO
509 ENDIF
510 ENDIF
511c
512C unloading with damage based on the energy.
513c
514 IF (iflag == 2) THEN
515 IF (nfunc == 1) THEN ! quasi-static only => 1D interpolation
516 yfac1 = uparam(nfunc + 11)
517 DO i=1,nel
518 ecurent(i)= zero
519 emax(i) = zero
520 emin(i) = ep20
521 et(i) = zero
522 ENDDO
523!
524 DO k=1,3
525 kk = (k-1)*nfunc + 1
526 DO i= 1,nindx
527 ii = indx(i)
528 iad(i) = npf(ifunc(1)) / 2 + 1
529 ipos(i) = vartmp(ii,kk)
530 ilen(i) = npf(ifunc(1)+1) / 2 - iad(i) - ipos(i)
531 ENDDO
532
533 CALL vinter2(tf,iad,ipos,ilen,nindx,strain(1,k),dydx,sqstat(1,k))
534
535 DO i=1,nindx
536 ii = indx(i)
537 sqstat(i,k) = yfac1 * sqstat(i,k)
538 IF (tflag == 2 .AND. strain(i,k) < zero) sqstat(i,k) = e0*strain(i,k)
539 s(i,k) = sqstat(i,k)
540 emax(i) = max(emax(i), yfac1*dydx(i))
541 emin(i) = min(emin(i), yfac1*dydx(i))
542C compute current energy
543 ecurent(i) = ecurent(i) + half*strain(i,k)*sqstat(i,k)
544 vartmp(ii,kk) = ipos(i)
545 END DO
546 END DO
547 DO i=1,nindx
548 ii = indx(i)
549 et(ii) = emin(i)/e0 ! not used
550 e(i) = max(e0,emax(i))
551 e_min(i) =emin(i)
552 yld(i) = sqrt(s(i,1)**2 + s(i,2)**2 + s(i,3)**2)
553 ENDDO
554c
555 ELSEIF (nfunc > 1) THEN ! strain rate dependent => 2D interpolation
556C
557 yfac1 = uparam(nfunc + 11)
558 DO i=1,nindx
559 ecurent(i) = zero
560 emax(i) = zero
561 emin(i) = ep20
562 END DO
563!
564 DO k=1,3
565 kk = (k-1)*nfunc + 1 ! static
566 DO i= 1,nindx
567 ii = indx(i)
568 iad(i) = npf(ifunc(1)) / 2 + 1
569 ipos(i) = vartmp(ii,kk)
570 ilen(i) = npf(ifunc(1)+1) / 2 - iad(i) - ipos(i)
571 ENDDO
572
573 CALL vinter2(tf,iad,ipos,ilen,nindx,strain(1,k),dydx,stmp1)
574
575 ! unloading + loading : 1:NE_UNL + NE_L
576 DO i=1,nindx
577 ii = indx(i)
578 sqstat(i,k) = yfac1*stmp1(i)
579 IF (tflag == 2 .AND. strain(i,k) < zero ) sqstat(i,k) = e0*strain(i,k)
580 s(i,k) = sqstat(i,k)
581 emax(i) = max(emax(i), yfac1*dydx(i))
582 emin(i) = min(emin(i), yfac1*dydx(i))
583 ! compute current energy
584 ecurent(i) = ecurent(i) + half*strain(i,k)*sqstat(i,k)
585 vartmp(ii,kk) = ipos(i)
586 END DO
587
588 jj(1:nindx) = 1
589 DO i=1,nindx
590 DO j = 2,nfunc - 1
591 IF (strainrate(i,k) >= uparam(10 + j)) jj(i) = j
592 ENDDO
593 END DO
594 DO i=1,nindx
595 rate(i,1) = uparam(10 + jj(i))
596 rate(i,2) = uparam(10 + jj(i)+1)
597 yfac(i,1) = uparam(10 + nfunc + jj(i))
598 yfac(i,2) = uparam(10 + nfunc + jj(i)+1)
599 END DO
600
601 DO i=1,nindx
602 ii = indx(i)
603 j1 = jj(i)
604 kk = (k-1)*nfunc + j1
605 ipos(i) = vartmp(ii,kk)
606 iad(i) = npf(ifunc(j1)) / 2 + 1
607 ilen(i) = npf(ifunc(j1)+1) / 2 - iad(i) - ipos(i)
608 END DO
609 CALL vinter2(tf,iad,ipos,ilen,nindx,strain(1,k),dydx,stmp1)
610 DO n=1,ne_l
611 i = indx_l(n)
612 ii = indx(i)
613 j1 = jj(i)
614 kk = (k-1)*nfunc + j1
615 vartmp(ii,kk) = ipos(i)
616 stmp1(i) = yfac(i,1) * stmp1(i)
617 emax(i) = max(emax(i), yfac(i,1)*dydx(i))
618 emin(i) = max(emin(i), yfac(i,1)*dydx(i))
619 END DO
620
621 DO i=1,nindx
622 ii = indx(i)
623 j2 = jj(i)+1
624 kk = (k-1)*nfunc + j2
625 ipos(i) = vartmp(ii,kk)
626 iad(i) = npf(ifunc(j2)) / 2 + 1
627 ilen(i) = npf(ifunc(j2)+1) / 2 - iad(i) - ipos(i)
628 END DO
629 CALL vinter2(tf,iad,ipos,ilen,nindx,strain(1,k),dydx,stmp2)
630 DO i=1,nindx
631 ii = indx(i)
632 j2 = jj(i) + 1
633 kk = (k-1)*nfunc + j2
634 vartmp(i,kk) = ipos(i)
635 stmp2(i) = yfac(i,2) * stmp2(i)
636 emax(i) = max(emax(i), yfac(i,2)*dydx(i))
637 emin(i) = max(emin(i), yfac(i,2)*dydx(i))
638 vartmp(ii,kk) = ipos(i)
639 END DO
640
641 DO ii=1,ne_l
642 i = indx_l(ii)
643 fac = (strainrate(i,k) - rate(i,1)) / (rate(i,2) - rate(i,1))
644 s(i,k) = stmp1(i) + fac*(stmp2(i) - stmp1(i))
645 IF (tflag == 2 .AND. strain(i,k) < zero ) s(i,k) = e0*strain(i,k)
646 END DO
647
648 ENDDO ! K
649!
650 DO i=1,nindx
651 ii = indx(i)
652 e(i) = max(e0,emax(i))
653 e(i) = emax(i)
654 et(ii) = emin(i)/e0 ! not used
655 yld(i) = sqrt(s(i,1)**2 + s(i,2)**2 + s(i,3)**2)
656 ENDDO
657
658 ENDIF ! NFUNC
659!-----------------------------------------------------------------------
660#include "vectorize.inc"
661 DO n=1,nindx
662 i = indx(n)
663 delta = epst(n) - uvar(i,6)
664 uvar(i,4) = uvar(i,4) + half*(yld(n) + uvar(i,1))*delta
665 uvar(i,4) = max(zero, uvar(i,4))
666 uvar(i,2) = max(uvar(i,2) , uvar(i,4))
667 uvar(i,1) = yld(n)
668 uvar(i,6) = epst(n)
669 ecurent(n) = uvar(i,4)
670 ENDDO
671C
672C idam is a hidden flag, only idam=0 is activated, Idam > 0 not tested.
673 IF (idam == 0) THEN
674#include "vectorize.inc"
675 DO n=1,ne_unl
676 ii = indx_unl(n)
677 i = indx(ii)
678 IF (uvar(i,2) > zero) THEN
679 dam = one - (ecurent(ii)/uvar(i,2))**shape
680 dam = dam**alpha
681 dam = one - (one - hys)*dam
682 uvar(i,7) = dam
683C global
684 DO k=1,3
685 s(ii,k)= dam*s(ii,k)
686 ENDDO
687 ENDIF
688 ENDDO ! NE_UNL
689c
690 ELSE
691
692C damage by direction to be tested for
693#include "vectorize.inc"
694 DO n=1,ne_unl
695 ii = indx_unl(n)
696 i = indx(ii)
697 IF (uvar(i,2) > zero) THEN
698 dam = one - (ecurent(ii)/uvar(i,2))**shape
699 dam = one - (one - hys)*dam
700 uvar(i,7) = dam
701 DO k=1,3
702 IF(iload(ii,k) < 0)s(ii,k) = dam*s(ii,k)
703 ENDDO
704 ENDIF
705 ENDDO ! nel
706 ENDIF ! IDAM
707 ENDIF ! iflag=2
708C
709C =====================================================
710 IF( fail > 0) THEN
711#include "vectorize.inc"
712 DO i = 1,nindx
713 ii = indx(i)
714C S > 0 for compression - curve definition
715C S < 0 for traction
716 t1 = -s(i,1)
717 t2 = -s(i,2)
718 t3 = -s(i,3)
719 IF(t1 >= tcut0 .OR. t2 >= tcut0 .OR. t3 >= tcut0 ) off(ii) = four_over_5
720 t1 = t1/ev(i,2)/ev(i,3)
721 t2 = t2/ev(i,1)/ev(i,3)
722 t3 = t3/ev(i,1)/ev(i,2)
723C
724C cauchy to glabale
725C
726 signxx(ii) = dirprv(i,1,1)*dirprv(i,1,1)*t1
727 . + dirprv(i,1,2)*dirprv(i,1,2)*t2
728 . + dirprv(i,1,3)*dirprv(i,1,3)*t3
729
730 signyy(ii) = dirprv(i,2,2)*dirprv(i,2,2)*t2
731 . + dirprv(i,2,3)*dirprv(i,2,3)*t3
732 . + dirprv(i,2,1)*dirprv(i,2,1)*t1
733
734 signzz(ii) = dirprv(i,3,3)*dirprv(i,3,3)*t3
735 . + dirprv(i,3,1)*dirprv(i,3,1)*t1
736 . + dirprv(i,3,2)*dirprv(i,3,2)*t2
737 signxy(ii) = dirprv(i,1,1)*dirprv(i,2,1)*t1
738 . + dirprv(i,1,2)*dirprv(i,2,2)*t2
739 . + dirprv(i,1,3)*dirprv(i,2,3)*t3
740
741 signyz(ii) = dirprv(i,2,2)*dirprv(i,3,2)*t2
742 . + dirprv(i,2,3)*dirprv(i,3,3)*t3
743 . + dirprv(i,2,1)*dirprv(i,3,1)*t1
744
745 signzx(ii) = dirprv(i,3,3)*dirprv(i,1,3)*t3
746 . + dirprv(i,3,1)*dirprv(i,1,1)*t1
747 . + dirprv(i,3,2)*dirprv(i,1,2)*t2
748 !
749 ! computing E
750 e_old = uvar(ii,8)
751 epss = epst(i) - yld(i)/ e_old
752 epss = max(zero, epss)
753 epss = min(one, epss)
754 aa = e_max - e0
755 de = aa*(epss - uvar(ii,10))
756 e_new = aa*epss + e0
757 e_new= min(e_max, e_new)
758 uvar(ii,10) = epss
759 uvar(ii,8) = e_new
760 aa = e_new*(one-nu)/(one + nu)/(one - two*nu)
761 soundsp(ii) = sqrt(aa/rho0(ii))
762 viscmax(ii) = zero
763 ENDDO ! NINDX
764 ELSEIF(fail == 0) THEN
765#include "vectorize.inc"
766 DO i = 1,nindx
767 ii = indx(i)
768C S > 0 for compression - curve definition
769C S < 0 for traction
770 t1 = -s(i,1)
771 t2 = -s(i,2)
772 t3 = -s(i,3)
773 tcut = tcut0
774 IF(unload(i) == 1 ) tcut = uvar(ii,7)*tcut0
775 t1 = min(tcut, t1)
776 t2 = min(tcut, t2)
777 t3 = min(tcut, t3)
778 t1 = t1/ev(i,2)/ev(i,3)
779 t2 = t2/ev(i,1)/ev(i,3)
780 t3 = t3/ev(i,1)/ev(i,2)
781C
782C cauchy to glabale
783C
784 signxx(ii) = dirprv(i,1,1)*dirprv(i,1,1)*t1
785 . + dirprv(i,1,2)*dirprv(i,1,2)*t2
786 . + dirprv(i,1,3)*dirprv(i,1,3)*t3
787
788 signyy(ii) = dirprv(i,2,2)*dirprv(i,2,2)*t2
789 . + dirprv(i,2,3)*dirprv(i,2,3)*t3
790 . + dirprv(i,2,1)*dirprv(i,2,1)*t1
791
792 signzz(ii) = dirprv(i,3,3)*dirprv(i,3,3)*t3
793 . + dirprv(i,3,1)*dirprv(i,3,1)*t1
794 . + dirprv(i,3,2)*dirprv(i,3,2)*t2
795 signxy(ii) = dirprv(i,1,1)*dirprv(i,2,1)*t1
796 . + dirprv(i,1,2)*dirprv(i,2,2)*t2
797 . + dirprv(i,1,3)*dirprv(i,2,3)*t3
798
799 signyz(ii) = dirprv(i,2,2)*dirprv(i,3,2)*t2
800 . + dirprv(i,2,3)*dirprv(i,3,3)*t3
801 . + dirprv(i,2,1)*dirprv(i,3,1)*t1
802
803 signzx(ii) = dirprv(i,3,3)*dirprv(i,1,3)*t3
804 . + dirprv(i,3,1)*dirprv(i,1,1)*t1
805 . + dirprv(i,3,2)*dirprv(i,1,2)*t2
806 ! computing E
807 e_old = uvar(ii,8)
808 epss = epst(i) - yld(i)/ e_old
809 epss = max(zero, epss)
810 epss = min(one, epss)
811 aa = e_max - e0
812 de = aa*(epss - uvar(ii,10))
813 e_new = aa*epss + e0
814 e_new= min(e_max, e_new)
815 uvar(ii,10) = epss
816 uvar(ii,8) = e_new
817 aa = e_new*(one-nu)/(one + nu)/(one - two*nu)
818 soundsp(ii) = sqrt(aa/rho0(ii))
819 viscmax(ii) = zero
820 ENDDO ! NINDX
821 ENDIF ! FAIL
822C==================================================
823 IF (impl_s > 0 .OR. ihet > 1) THEN
824#include "vectorize.inc"
825 DO n=1,nindx
826 i = indx(n)
827 et(i) = yld(n)/max(em20,epst(n))
828 et(i) = min(one , et(i)/e_max)
829 IF(et(i) == zero) et(i) = one
830 ENDDO
831 ENDIF
832c
833 1000 FORMAT(1x,'TCUT REACHED, DELETED SOLID ELEMENT ',i10)
834 1100 FORMAT(1x,'TCUT REACHED, DELETED SOLID ELEMENT ',i10,1x,'AT TIME :',1pe12.4)
835c
836c-----------
837C------------------------------------
838 RETURN
839 END SUBROUTINE sigeps90
840C
sigeps90
subroutine sigeps90(nel, nuvar, nfunc, ifunc, npf, tf, time, uparam, rho0, epspxx, epspyy, epspzz, epspxy, epspyz, epspzx, epsxx, epsyy, epszz, epsxy, epsyz, epszx, sigoxx, sigoyy, sigozz, sigoxy, sigoyz, sigozx, signxx, signyy, signzz, signxy, signyz, signzx, soundsp, viscmax, uvar, nvartmp, vartmp, ismstr, israte, asrate, offg, ihet, et, epsd, off, ngl)
Definition sigeps90.F:44
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
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
vinter2
subroutine vinter2(tf, iad, ipos, ilen, nel0, x, dydx, y)
Definition vinter.F:143