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 ../starter/source/materials/mat/mat090/sigeps90.F
25!||--- called by ------------------------------------------------------
26!|| mulaw ../starter/source/materials/mat_share/mulaw.F
27!||--- calls -----------------------------------------------------
28!|| finter ../starter/source/tools/curve/finter.F
29!|| valpvec ../starter/source/materials/tools/matrix.F
30!|| valpvecdp ../starter/source/materials/tools/matrix.F
31!||====================================================================
32 SUBROUTINE sigeps90(
33 1 NEL ,NUPARAM,NUVAR ,NFUNC ,IFUNC ,NPF ,
34 2 TF ,TIME ,TIMESTEP,UPARAM ,RHO0 ,RHO ,
35 3 VOLUME ,EINT ,
36 4 EPSPXX ,EPSPYY ,EPSPZZ ,EPSPXY ,EPSPYZ ,EPSPZX ,
37 5 DEPSXX ,DEPSYY ,DEPSZZ ,DEPSXY ,DEPSYZ ,DEPSZX ,
38 6 EPSXX ,EPSYY ,EPSZZ ,EPSXY ,EPSYZ ,EPSZX ,
39 7 SIG0XX ,SIG0YY ,SIG0ZZ ,SIG0XY ,SIG0YZ ,SIG0ZX ,
40 8 SIGNXX ,SIGNYY ,SIGNZZ ,SIGNXY ,SIGNYZ ,SIGNZX ,
41 A SOUNDSP,VISCMAX,UVAR ,OFF ,NGL ,ISMSTR )
42C-----------------------------------------------
43C I m p l i c i t T y p e s
44C-----------------------------------------------
45#include "implicit_f.inc"
46C-----------------------------------------------
47C G l o b a l P a r a m e t e r s
48C-----------------------------------------------
49C---------+---------+---+---+--------------------------------------------
50C VAR | SIZE |TYP| RW| DEFINITION
51C---------+---------+---+---+--------------------------------------------
52C NEL | 1 | I | R | SIZE OF THE ELEMENT GROUP NEL
53C NUPARAM | 1 | I | R | SIZE OF THE USER PARAMETER ARRAY
54C NUVAR | 1 | I | R | NUMBER OF USER ELEMENT VARIABLES
55C---------+---------+---+---+--------------------------------------------
56C MFUNC | 1 | I | R | NUMBER FUNCTION USED FOR THIS USER LAW not used
57C KFUNC | NFUNC | I | R | FUNCTION INDEX not used
58C NPF | * | I | R | FUNCTION ARRAY
59C TF | * | F | R | FUNCTION ARRAY
60C---------+---------+---+---+--------------------------------------------
61C TIME | 1 | F | R | CURRENT TIME
62C TIMESTEP| 1 | F | R | CURRENT TIME STEP
63C UPARAM | NUPARAM | F | R | USER MATERIAL PARAMETER ARRAY
64C RHO0 | NEL | F | R | INITIAL DENSITY
65C RHO | NEL | F | R | DENSITY
66C VOLUME | NEL | F | R | VOLUME
67C EINT | NEL | F | R | TOTAL INTERNAL ENERGY
68C EPSPXX | NEL | F | R | STRAIN RATE XX
69C EPSPYY | NEL | F | R | STRAIN RATE YY
70C ... | | | |
71C DEPSXX | NEL | F | R | STRAIN INCREMENT XX
72C DEPSYY | NEL | F | R | STRAIN INCREMENT YY
73C ... | | | |
74C EPSXX | NEL | F | R | STRAIN XX
75C EPSYY | NEL | F | R | STRAIN YY
76C ... | | | |
77C SIG0XX | NEL | F | R | OLD ELASTO PLASTIC STRESS XX
78C SIG0YY | NEL | F | R | OLD ELASTO PLASTIC STRESS 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"
94C
95 INTEGER, INTENT(IN) :: NEL, NUPARAM, NUVAR, NGL(NEL),ISMSTR
96
97 my_real, INTENT(IN) ::
98 . TIME,TIMESTEP,UPARAM(*),
99 . RHO(NEL),RHO0(NEL),VOLUME(NEL),EINT(NEL),
100 . EPSPXX(NEL),EPSPYY(NEL),EPSPZZ(NEL),
101 . EPSPXY(NEL),EPSPYZ(NEL),EPSPZX(NEL),
102 . DEPSXX(NEL),DEPSYY(NEL),DEPSZZ(NEL),
103 . DEPSXY(NEL),DEPSYZ(NEL),DEPSZX(NEL),
104 . epsxx(nel) ,epsyy(nel) ,epszz(nel) ,
105 . epsxy(nel) ,epsyz(nel) ,epszx(nel) ,
106 . sig0xx(nel),sig0yy(nel),sig0zz(nel),
107 . sig0xy(nel),sig0yz(nel),sig0zx(nel)
108C-----------------------------------------------
109C O U T P U T A r g u m e n t s
110C-----------------------------------------------
111 my_real,INTENT(OUT)::
112 . signxx(nel),signyy(nel),signzz(nel),
113 . signxy(nel),signyz(nel),signzx(nel),
114 . soundsp(nel),viscmax(nel)
115C-----------------------------------------------
116C I N P U T O U T P U T A r g u m e n t s
117C-----------------------------------------------
118 my_real,INTENT(INOUT) ::
119 . uvar(nel,nuvar), off(nel)
120C-----------------------------------------------
121C VARIABLES FOR FUNCTION INTERPOLATION
122C-----------------------------------------------
123 INTEGER, INTENT(IN) :: NPF(*), NFUNC, IFUNC(NFUNC)
124 my_real, INTENT(IN) :: TF(*)
125 my_real, EXTERNAL :: FINTER
126C EXTERNAL FINTER
127C Y = FINTER(IFUNC(J),X,NPF,TF,DYDX)
128C Y : y = f(x)
129C X : x
130C DYDX : f'(x) = dy/dx
131C IFUNC(J): FUNCTION INDEX
132C J : FIRST(J=1), SECOND(J=2) .. FUNCTION USED FOR THIS LAW
133C NPF,TF : FUNCTION PARAMETER
134C-----------------------------------------------
135C L o c a l V a r i a b l e s
136C-----------------------------------------------
137 INTEGER
138 . II,I,J,K,I1,J1,J2,IFLAG,ILOAD(NEL,3),ILOADE(NEL),
139 . IDAM,INDX_L(NEL),INDX_UNL(NEL),NE_L,NE_UNL
140 my_real
141 . E0,AA,EPSMAX,G,NU,SHAPE,HYS,EMAX,
142 . YFAC,YFACJ1,YFACJ2,RATEJ1,RATEJ2, EPSE,EP1,
143 . EP2,EP3,EP4,EP5,EP6,ERT11,ERT12,ERT13,ERT21,
144 . ERT22,ERT23,ERT31,ERT32,ERT33,SJ1,SJ2,FAC,T1,T2,T3,
145 . DAM,EPE,EMIN,E_MIN(NEL),E1,E2,TMP,QUASI_EINT
146 my_real
147 . av(6,nel), evv(3,nel),ev(nel,3),strain(nel,3),
148 . strainrate(nel,3),s(nel,3),sqstat(nel,3),df(3),
149 . dirprv(3,3,nel),epsp(3),ecurent(nel),e(nel),deint,
150 . rateeps
151C-----------------------------------------------
152C USER VARIABLES INITIALIZATION
153C-----------------------------------------------
154 e0 = uparam(1)
155 g = uparam(4)
156 nu = uparam(5)
157 shape = uparam(6)
158 hys = uparam(7)
159 iflag = uparam(9)
160 idam = uparam(10)
161C
162C-----------------------------------------------
163C
164 DO i=1,nel
165 av(1,i) = epsxx(i)
166 av(2,i) = epsyy(i)
167 av(3,i) = epszz(i)
168 av(4,i) = half*epsxy(i)
169 av(5,i) = half*epsyz(i)
170 av(6,i) = half*epszx(i)
171 ENDDO
172CEigenvalues needed to be calculated in double precision
173C for a simple precision executing*
174 IF (iresp==1) THEN
175 CALL valpvecdp(av,evv,dirprv,nel)
176 ELSE
177 CALL valpvec(av,evv,dirprv,nel)
178 ENDIF
179C-ISMSTR=0-NO SMALL STRAIN OPTION:STRAINS ARE LOGARITHMIC, STRESS IS CAUCHY
180C-ISMSTR=1-SMALL STRAIN OPTION:STRAINS ARE ENGINEERING, STRESS IS CAUCHY
181C-ISMSTR=2-SMALL STRAIN OPTION:STRAINS ARE ENGINEERING, STRESS IS BIOT
182C-ISMSTR=3-NO SMALL STRAIN OPTION:STRESS IS BIOT
183 IF(ismstr==0.OR.ismstr==2.OR.ismstr==4) THEN
184 DO i=1,nel
185C ---- (STRAIN IS LOGARITHMIC)
186 ev(i,1)=exp(evv(1,i))
187 ev(i,2)=exp(evv(2,i))
188 ev(i,3)=exp(evv(3,i))
189 ENDDO
190 ELSEIF(ismstr==10.OR.ismstr==12) THEN
191 DO i =1,nel
192 ev(i,1)=sqrt(evv(1,i) + one )
193 ev(i,2)=sqrt(evv(2,i) + one )
194 ev(i,3)=sqrt(evv(3,i) + one )
195 ENDDO
196 ELSE
197C ---- STRAIN IS ENGINEERING)
198 DO i=1,nel
199 ev(i,1)=evv(1,i) + one
200 ev(i,2)=evv(2,i) + one
201 ev(i,3)=evv(3,i) + one
202 ENDDO
203 ENDIF
204
205C engineering strain and strain rate
206 DO i=1,nel
207C engineering strain e = lambda-1 , according the input curve
208C e=1-lambda (e > 0 compression and e < 0 traction)
209 strain(i,1) = one - ev(i,1)
210 strain(i,2) = one - ev(i,2)
211 strain(i,3) = one - ev(i,3)
212C
213 ep1 = epspxx(i)
214 ep2 = epspyy(i)
215 ep3 = epspzz(i)
216 ep4 = half*epspxy(i)
217 ep5 = half*epspyz(i)
218 ep6 = half*epspzx(i)
219C phi_trans*L*phi_t
220 ert11 =dirprv(1,1,i)*ep1 + dirprv(2,1,i)*ep4 + dirprv(3,1,i)*ep6
221 ert12 =dirprv(1,2,i)*ep1 + dirprv(2,2,i)*ep4 + dirprv(3,2,i)*ep6
222 ert13 =dirprv(1,3,i)*ep1 + dirprv(2,3,i)*ep4 + dirprv(3,3,i)*ep6
223
224 ert21 =dirprv(1,1,i)*ep4 + dirprv(2,1,i)*ep2 + dirprv(3,1,i)*ep5
225 ert22 =dirprv(1,2,i)*ep4 + dirprv(2,2,i)*ep2 + dirprv(3,2,i)*ep5
226 ert23 =dirprv(1,3,i)*ep4 + dirprv(2,3,i)*ep2 + dirprv(3,3,i)*ep5
227
228 ert31 =dirprv(1,1,i)*ep6 + dirprv(2,1,i)*ep5 + dirprv(3,1,i)*ep3
229 ert32 =dirprv(1,2,i)*ep6 + dirprv(2,2,i)*ep5 + dirprv(3,2,i)*ep3
230 ert33 =dirprv(1,3,i)*ep6 + dirprv(2,3,i)*ep5 + dirprv(3,3,i)*ep3
231C
232 epsp(1) = dirprv(1,1,i)*ert11 + dirprv(2,1,i)*ert21
233 . + dirprv(3,1,i)*ert31
234 epsp(2) = dirprv(1,2,i)*ert12 + dirprv(2,2,i)*ert22
235 . + dirprv(3,2,i)*ert32
236 epsp(3) = dirprv(1,3,i)*ert13 + dirprv(2,3,i)*ert23
237 . + dirprv(3,3,i)*ert33
238 strainrate(i,1) = abs(epsp(1))
239 strainrate(i,2) = abs(epsp(2))
240 strainrate(i,3) = abs(epsp(3))
241C computing energy increase
242!! YFAC = UPARAM(NFUNC + 11)
243!! QUASI_EINT= ZERO
244!! DO K=1,3
245!! EPSE = STRAIN(I,K)
246!! SQSTAT(I,K) = YFAC*FINTER(IFUNC(1),EPSE,NPF,TF,DF(K))
247C compute current energy
248!! QUASI_EINT= QUASI_EINT + HALF*STRAIN(I,K)*SQSTAT(I,K)
249!! ENDDO
250!! DEINT = QUASI_EINT - UVAR(I,9)
251!! UVAR(I,9) = QUASI_EINT
252C -----------
253C check loading and unloading.
254 iloade(i) = 1
255 DO k=1,3
256 epe = epsp(k)*strain(i,k)
257 iload(i,k) = 1
258 IF(epe > em10 )iload(i,k) = -1
259 IF(iload(i,k) == -1 ) iloade(i) = -1
260 ENDDO
261C
262 DO k=1,3
263 IF(iload(i,k) == -1) strainrate(i,k) = zero
264 ENDDO
265
266 uvar(i,3) = strainrate(i,1)
267 uvar(i,4) = strainrate(i,2)
268 uvar(i,5) = strainrate(i,3)
269C
270 rateeps = max(strainrate(i,1),strainrate(i,2),strainrate(i,3))
271
272!! EPSD(I) = RATEEPS
273 ENDDO
274C sous groupe
275 indx_l(1:nel) = 0
276 indx_unl(1:nel) = 0
277 ne_l = 0
278 ne_unl = 0
279C
280 DO i=1,nel
281 IF(iloade(i) == 1) THEN
282 ne_l = ne_l + 1
283 indx_l(ne_l) = i
284 uvar(i,8) = one
285 ELSEIF(iloade(i) == -1 ) THEN
286 ne_unl = ne_unl + 1
287 indx_unl(ne_unl) = i
288 uvar(i,8) = -one
289 ENDIF
290 ENDDO
291C case with unloading stress-strain curve only the quasistatic curve is used
292 IF(iflag == 1) THEN
293 DO i=1,nel
294 yfac = uparam(nfunc + 11)
295 emax = zero
296 DO k=1,3
297 epse = strain(i,k)
298 s(i,k) = yfac*finter(ifunc(1),epse,npf,tf,df(k))
299 emax = max(emax, yfac*df(k))
300 ENDDO
301 e(i) = max(e0,emax)
302 ENDDO
303 ENDIF
304C unloading with damage based on the energy.
305 IF(iflag == 2) THEN
306 DO i=1,nel
307 yfac = uparam(nfunc + 11)
308 ecurent(i)= zero
309 emax = zero
310 emin = ep20
311 DO k=1,3
312 epse = strain(i,k)
313 sqstat(i,k) = yfac*finter(ifunc(1),epse,npf,tf,df(k))
314 s(i,k) = sqstat(i,k)
315 emax = max(emax, yfac*df(k))
316 emin = min(emin, yfac*df(k))
317C computecurent energy
318 ecurent(i)= ecurent(i) + half*strain(i,k)*sqstat(i,k)
319!! ECURENT(I)= ECURENT(I) + HALF*EV(I,K)*ABS(SQSTAT(I,K))
320 ENDDO
321 uvar(i,2) = max(uvar(i,2),ecurent(i))
322 e(i) = max(e0,emax)
323 e_min(i) =emin
324 ENDDO
325C flag is a hidden flag, only idam=0 is activated, Idam > o not tested.
326 IF(idam == 0) THEN
327 DO ii=1,ne_unl
328 i = indx_unl(ii)
329 uvar(i,1) = ecurent(i)
330 IF(uvar(i,2) > zero) THEN
331 dam = one - (ecurent(i)/uvar(i,2))**shape
332 dam = one - (one - hys)*dam
333C global
334 DO k=1,3
335 s(i,k)= dam*s(i,k)
336 ENDDO
337 ENDIF
338 ENDDO ! NE_UNL
339 ELSE ! not tested Idam/=0
340C damage by direction to be tested for
341 DO ii=1,ne_unl
342 i = indx_unl(ii)
343 uvar(i,1) = ecurent(i)
344 IF(uvar(i,2) > zero) THEN
345 dam = one - (ecurent(i)/uvar(i,2))**shape
346 dam = one - (one - hys)*dam
347 DO k=1,3
348 IF(iload(i,k) < 0)s(i,k) = dam*s(i,k)
349 ENDDO
350 ENDIF
351 ENDDO ! nel
352 ENDIF ! IDAM
353 ENDIF ! iflag=2
354C
355C =====================================================
356 DO i = 1,nel
357C S > 0 for compression - curve definition
358C S < 0 for traction
359 t1 = -s(i,1)/ev(i,2)/ev(i,3)
360 t2 = -s(i,2)/ev(i,1)/ev(i,3)
361 t3 = -s(i,3)/ev(i,1)/ev(i,2)
362C
363C cauchy to glabale
364C
365 signxx(i) = dirprv(1,1,i)*dirprv(1,1,i)*t1
366 . + dirprv(1,2,i)*dirprv(1,2,i)*t2
367 . + dirprv(1,3,i)*dirprv(1,3,i)*t3
368
369 signyy(i) = dirprv(2,2,i)*dirprv(2,2,i)*t2
370 . + dirprv(2,3,i)*dirprv(2,3,i)*t3
371 . + dirprv(2,1,i)*dirprv(2,1,i)*t1
372
373 signzz(i) = dirprv(3,3,i)*dirprv(3,3,i)*t3
374 . + dirprv(3,1,i)*dirprv(3,1,i)*t1
375 . + dirprv(3,2,i)*dirprv(3,2,i)*t2
376
377 signxy(i) = dirprv(1,1,i)*dirprv(2,1,i)*t1
378 . + dirprv(1,2,i)*dirprv(2,2,i)*t2
379 . + dirprv(1,3,i)*dirprv(2,3,i)*t3
380
381 signyz(i) = dirprv(2,2,i)*dirprv(3,2,i)*t2
382 . + dirprv(2,3,i)*dirprv(3,3,i)*t3
383 . + dirprv(2,1,i)*dirprv(3,1,i)*t1
384
385 signzx(i) = dirprv(3,3,i)*dirprv(1,3,i)*t3
386 . + dirprv(3,1,i)*dirprv(1,1,i)*t1
387 . + dirprv(3,2,i)*dirprv(1,2,i)*t2
388C==================================================
389 aa = e(i)*(one-nu)/(one + nu)/(one - two*nu)
390 soundsp(i) = sqrt(aa/rho0(i))
391 viscmax(i) = zero
392C
393 ENDDO
394CE0
395C------------------------------------
396 RETURN
397 END
398C
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
valpvecdp
subroutine valpvecdp(sig, val, vec, nel)
Definition matrix.F:533
valpvec
subroutine valpvec(sig, val, vec, nel)
Definition matrix.F:215
sigeps90
subroutine sigeps90(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, sig0xx, sig0yy, sig0zz, sig0xy, sig0yz, sig0zx, signxx, signyy, signzz, signxy, signyz, signzx, soundsp, viscmax, uvar, off, ngl, ismstr)
Definition sigeps90.F:42