OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps64c.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!|| sigeps64c ../engine/source/materials/mat/mat064/sigeps64c.F
25!||--- called by ------------------------------------------------------
26!|| mulawc ../engine/source/materials/mat_share/mulawc.F90
27!||--- calls -----------------------------------------------------
28!|| finter ../engine/source/tools/curve/finter.F
29!|| vinter ../engine/source/tools/curve/vinter.F
30!||====================================================================
31 SUBROUTINE sigeps64c(
32 1 NEL0, NUPARAM, NUVAR, MFUNC,
33 2 KFUNC, NPF, NPT0, IPT,
34 3 IFLAG, TF, TIME, TIMESTEP,
35 4 UPARAM, RHO0, AREA, EINT,
36 5 THKLY, EPSPXX, EPSPYY, EPSPXY,
37 6 EPSPYZ, EPSPZX, DEPSXX, DEPSYY,
38 7 DEPSXY, DEPSYZ, DEPSZX, EPSXX,
39 8 EPSYY, EPSXY, EPSYZ, EPSZX,
40 9 SIGOXX, SIGOYY, SIGOXY, SIGOYZ,
41 A SIGOZX, SIGNXX, SIGNYY, SIGNXY,
42 B SIGNYZ, SIGNZX, SIGVXX, SIGVYY,
43 C SIGVXY, SIGVYZ, SIGVZX, SOUNDSP,
44 D VISCMAX, THK, PLA, UVAR,
45 E OFF, NGL, IPM, MAT,
46 F ETSE, GS, VOL, YLD,
47 G TEMPEL, DIE, INLOC,
48 H DPLANL, JTHE, LOFF)
49C-----------------------------------------------
50C I m p l i c i t T y p e s
51C-----------------------------------------------
52#include "implicit_f.inc"
53C-----------------------------------------------
54C G l o b a l P a r a m e t e r s
55C-----------------------------------------------
56#include "mvsiz_p.inc"
57C---------+---------+---+---+--------------------------------------------
58C VAR | SIZE |TYP| RW| DEFINITION
59C---------+---------+---+---+--------------------------------------------
60C NEL0 | 1 | I | R | SIZE OF THE ELEMENT GROUP NEL0 F
61C NUPARAM | 1 | I | R | SIZE OF THE USER PARAMETER ARRAY
62C NUVAR | 1 | I | R | NUMBER OF USER ELEMENT VARIABLES
63C---------+---------+---+---+--------------------------------------------
64C NFUNC | 1 | I | R | NUMBER FUNCTION USED FOR THIS USER LAW
65C IFUNC | NFUNC | I | R | FUNCTION INDEX
66C NPF | * | I | R | FUNCTION ARRAY
67C NPT0 | 1 | I | R | NUMBER OF LAYERS OR INTEGRATION POINTS
68C IPT | 1 | I | R | LAYER OR INTEGRATION POINT NUMBER
69C IFLAG | * | I | R | GEOMETRICAL FLAGS
70C TF | * | F | R | FUNCTION ARRAY
71C---------+---------+---+---+--------------------------------------------
72C TIME | 1 | F | R | CURRENT TIME
73C TIMESTEP| 1 | F | R | CURRENT TIME STEP
74C UPARAM | NUPARAM | F | R | USER MATERIAL PARAMETER ARRAY
75C RHO0 | NEL0 | F | R | INITIAL DENSITY
76C AREA | NEL0 | F | R | AREA
77C EINT | 2*NEL0 | F | R | INTERNAL ENERGY(MEMBRANE,BENDING)
78C THKLY | NEL0 | F | R | LAYER THICKNESS
79C EPSPXX | NEL0 | F | R | STRAIN RATE XX
80C EPSPYY | NEL0 | F | R | STRAIN RATE YY
81C ... | | | |
82C DEPSXX | NEL0 | F | R | STRAIN INCREMENT XX
83C DEPSYY | NEL0 | F | R | STRAIN INCREMENT YY
84C ... | | | |
85C EPSXX | NEL0 | F | R | STRAIN XX
86C EPSYY | NEL0 | F | R | STRAIN YY
87C ... | | | |
88C SIGOXX | NEL0 | F | R | OLD ELASTO PLASTIC STRESS XX
89C SIGOYY | NEL0 | F | R | OLD ELASTO PLASTIC STRESS YY
90C ... | | | |
91C---------+---------+---+---+--------------------------------------------
92C SIGNXX | NEL0 | F | W | NEW ELASTO PLASTIC STRESS XX
93C SIGNYY | NEL0 | F | W | NEW ELASTO PLASTIC STRESS YY
94C ... | | | |
95C SIGVXX | NEL0 | F | W | VISCOUS STRESS XX
96C SIGVYY | NEL0 | F | W | VISCOUS STRESS YY
97C ... | | | |
98C SOUNDSP | NEL0 | F | W | SOUND SPEED (NEEDED FOR TIME STEP)
99C VISCMAX | NEL0 | F | W | MAXIMUM DAMPING MODULUS(NEEDED FOR TIME STEP)
100C---------+---------+---+---+--------------------------------------------
101C THK | NEL0 | F |R/W| THICKNESS
102C PLA | NEL0 | F |R/W| PLASTIC STRAIN
103C UVAR |NEL0*NUVAR| F |R/W| USER ELEMENT VARIABLE ARRAY
104C OFF | NEL0 | F |R/W| DELETED ELEMENT FLAG (=1. ON, =0. OFF)
105C---------+---------+---+---+--------------------------------------------
106C C o m m o n B l o c k s
107C-----------------------------------------------
108#include "param_c.inc"
109#include "com01_c.inc"
110C-----------------------------------------------
111C I N P U T A r g u m e n t s
112C-----------------------------------------------
113 INTEGER, INTENT(IN) :: JTHE
114 INTEGER NEL0, NUPARAM, NUVAR, NPT0, IPT,IFLAG(*),
115 . NGL(NEL0),MAT(NEL0),IPM(NPROPMI,*),INLOC
116 my_real TIME,TIMESTEP,UPARAM(*),
117 . AREA(NEL0),RHO0(NEL0),EINT(NEL0,2),
118 . THKLY(NEL0),PLA(NEL0),
119 . EPSPXX(NEL0),EPSPYY(NEL0),
120 . EPSPXY(NEL0),EPSPYZ(NEL0),EPSPZX(NEL0),
121 . DEPSXX(NEL0),DEPSYY(NEL0),
122 . DEPSXY(NEL0),DEPSYZ(NEL0),DEPSZX(NEL0),
123 . EPSXX(NEL0) ,EPSYY(NEL0) ,
124 . EPSXY(NEL0) ,EPSYZ(NEL0) ,EPSZX(NEL0) ,
125 . SIGOXX(NEL0),SIGOYY(NEL0),
126 . SIGOXY(NEL0),SIGOYZ(NEL0),SIGOZX(NEL0),
127 . GS(*),VOL(NEL0),TEMPEL(NEL0),
128 . DIE(NEL0),DPLANL(NEL0)
129 my_real, DIMENSION(NEL0), INTENT(IN) :: LOFF
130C-----------------------------------------------
131C O U T P U T A r g u m e n t s
132C-----------------------------------------------
133 my_real
134 . signxx(nel0),signyy(nel0),
135 . signxy(nel0),signyz(nel0),signzx(nel0),
136 . sigvxx(nel0),sigvyy(nel0),
137 . sigvxy(nel0),sigvyz(nel0),sigvzx(nel0),
138 . soundsp(nel0),viscmax(nel0),etse(nel0)
139C-----------------------------------------------
140C I N P U T O U T P U T A r g u m e n t s
141C-----------------------------------------------
142 my_real uvar(nel0,nuvar), off(nel0),thk(nel0),yld(nel0)
143C-----------------------------------------------
144C VARIABLES FOR FUNCTION INTERPOLATION
145C-----------------------------------------------
146 INTEGER NPF(*), MFUNC, KFUNC(MFUNC)
147 my_real FINTER ,TF(*)
148 EXTERNAL FINTER
149C Y = FINTER(IFUNC(J),X,NPF,TF,DYDX)
150C Y : y = f(x)
151C X : x
152C DYDX : f'(x) = dy/dx
153C IFUNC(J): FUNCTION INDEX
154C J : FIRST(J=1), SECOND(J=2) .. FUNCTION USED FOR THIS LAW
155C NPF,TF : FUNCTION PARAMETER
156C-----------------------------------------------
157C L o c a l V a r i a b l e s
158C-----------------------------------------------
159 INTEGER I,J,N,NINDX,
160 . NMAX,INDEX(MVSIZ),
161 . IAD1(MVSIZ),IPOS1(MVSIZ),ILEN1(MVSIZ),
162 . IAD2(MVSIZ),IPOS2(MVSIZ),ILEN2(MVSIZ),
163 . nfunc, ifunc(2),mx
164 my_real
165 . e,a1,a2,g,g3,
166 . aa(mvsiz),bb(mvsiz),pp(mvsiz),qq(mvsiz),h(mvsiz),
167 .
168 . cn,
169 . v0, temp(mvsiz),vmc,fac,
170 . vol0,dvm, vm(mvsiz),nu,
171 . nnu2,nu1,nu2,nu3,nu4,
172 . nu5,nu6,svm(mvsiz),
173 . dydx1(mvsiz),dydx2(mvsiz), y1(mvsiz),y2(mvsiz)
174 my_real
175 . dpla(mvsiz),umr, r, cp,dezz,s1, s2,
176 . s3, dpla_j(mvsiz), yld_i,dr(mvsiz), p2,q2,nnu1,
177 . s11,s22,
178 . s12, vm2, a, b, c, sigz, f, df, p,md,
179 . cd,yfac(2),
180 . hl
181 DATA nmax/3/
182C-----------------------------------------------
183C USER VARIABLES INITIALIZATION
184C-----------------------------------------------
185C
186 mx = mat(1)
187 nfunc = ipm(10,mx)
188 DO j=1,nfunc
189 ifunc(j)=ipm(10+j,mx)
190 ENDDO
191 e = uparam(1)
192 a1 = uparam(2)
193 a2 = uparam(3)
194 g = uparam(4)
195 nu = uparam(5)
196 cd = uparam(6)
197 cn = uparam(7)
198 md = uparam(8)
199 v0 = uparam(9)
200 vmc = uparam(10)
201 cp = uparam(11)
202 yfac(1) = uparam(13)
203 yfac(2) = uparam(14)
204 hl = uparam(15)
205 g3 = three*g
206 nnu1 = nu / (one - nu)
207 nnu2 = nnu1*nnu1
208 nu1 = one/(one - nu)
209 nu2 = one/(one + nu)
210 nu3 = one -nnu1
211 nu4 = one + nnu2 + nnu1
212 nu5 = one + nnu2 - two*nnu1
213 nu6 = half - nnu2 + half*nnu1
214 DO i=1,nel0
215 temp(i) = uparam(12)
216 ENDDO
217C-----------------------------------------------
218C Calculation of temperature. (conduction or adiabatic)
219C--------------------
220 IF(jthe > 0 ) THEN
221 DO i=1,nel0
222 temp(i) = tempel(i)
223 ENDDO
224 ELSE
225 DO i=1,nel0
226 vol0 = vol(i) * rho0(i)
227 vm(i) = uvar(i,2)
228 temp(i) = temp(i) + (eint(i,1)+ eint(i,2) + vm(i)*hl) * cp / vol0
229 ENDDO
230 ENDIF
231C
232 DO i=1,nel0
233C
234 signxx(i)=sigoxx(i) + a1*depsxx(i) + a2*depsyy(i)
235 signyy(i)=sigoyy(i) + a2*depsxx(i) + a1*depsyy(i)
236 signxy(i)=sigoxy(i) + g *depsxy(i)
237 signyz(i)=sigoyz(i) + gs(i)*depsyz(i)
238 signzx(i)=sigozx(i) + gs(i)*depszx(i)
239 vm(i) = uvar(i,2)
240 dpla(i) = zero
241C ...
242
243 soundsp(i) = sqrt(a1/rho0(i))
244 viscmax(i) = zero
245 etse(i) = one
246C
247 ENDDO
248C compute YLD stress
249 DO i=1,nel0
250 ipos1(i) = nint(uvar(i,3))
251 iad1(i) = npf(ifunc(1)) / 2 + 1
252 ilen1(i) = npf(ifunc(1)+1) / 2 - iad1(i) - ipos1(i)
253 ipos2(i) = nint(uvar(i,4))
254 iad2(i) = npf(ifunc(2)) / 2 + 1
255 ilen2(i) = npf(ifunc(2)+1) / 2 - iad2(i) - ipos2(i)
256 ENDDO
257C
258 CALL vinter(tf,iad1,ipos1,ilen1,nel0,pla,dydx1,y1)
259 CALL vinter(tf,iad2,ipos2,ilen2,nel0,pla,dydx2,y2)
260C Yld stress
261 DO i=1,nel0
262 y1(i)=y1(i)*yfac(1)
263 y2(i)=y2(i)*yfac(2)
264 fac = vm(i)/vmc
265 yld(i) = y1(i) + fac*(y2(i)-y1(i))
266 yld(i) = max(yld(i),em20)
267C collapse
268 dydx1(i)=dydx1(i)*yfac(1)
269 dydx2(i)=dydx2(i)*yfac(2)
270 h(i) = (dydx1(i) + fac*(dydx2(i)-dydx1(i)))
271 uvar(i,3) = ipos1(i)
272 uvar(i,4) = ipos2(i)
273 ENDDO
274C-------------------
275C PROJECTION
276C-------------------
277 IF(iflag(1)==0)THEN
278C radial projection
279 DO i=1,nel0
280 svm(i)=sqrt(signxx(i)*signxx(i)
281 . +signyy(i)*signyy(i)
282 . -signxx(i)*signyy(i)
283 . +three*signxy(i)*signxy(i))
284 r = min(one,yld(i)/max(em20,svm(i)))
285 signxx(i)=signxx(i)*r
286 signyy(i)=signyy(i)*r
287 signxy(i)=signxy(i)*r
288 umr = one - r
289cc DPLA(I) = OFF(I)*SVM(I)*UMR/(G3(I)+H(I))
290 dpla(i) = off(i)*svm(i)*umr/e
291 pla(i) = pla(i) + dpla(i)
292 s1=half*(signxx(i)+signyy(i))
293 IF (inloc == 0) THEN
294 dezz = dpla(i) * half*(signxx(i)+signyy(i)) / yld(i)
295 dezz=-(depsxx(i)+depsyy(i))*nnu1-nu3*dezz
296 thk(i) = thk(i) + dezz*thkly(i)*off(i)
297 ENDIF
298 IF(r<one) etse(i)= h(i)/(h(i)+e)
299 ENDDO
300 ELSEIF(iflag(1)==1)THEN
301C-------------------------
302C VON MISES CRITERION
303C-------------------------
304 DO i=1,nel0
305 h(i) = max(zero,h(i))
306 s1=signxx(i)+signyy(i)
307 s2=signxx(i)-signyy(i)
308 s3=signxy(i)
309 aa(i)=fourth*s1*s1
310 bb(i)=three_over_4*s2*s2+three*s3*s3
311 svm(i)=sqrt(aa(i)+bb(i))
312 IF (inloc == 0) THEN
313 dezz = -(depsxx(i)+depsyy(i))*nnu1
314 thk(i) = thk(i) + dezz*thkly(i)*off(i)
315 ENDIF
316 ENDDO
317C-------------------------
318C GATHER PLASTIC FLOW
319C-------------------------
320 nindx=0
321 DO i=1,nel0
322 IF(svm(i)>yld(i).AND.off(i)==one) THEN
323 nindx=nindx+1
324 index(nindx)=i
325 ENDIF
326 ENDDO
327C
328 IF(nindx/=0) THEN
329C---------------------------
330C DEP EN CONTRAINTE PLANE
331C---------------------------
332 DO j=1,nindx
333 i=index(j)
334 dpla_j(i)=(svm(i)-yld(i))/(g3+h(i))
335 etse(i)= h(i)/(h(i)+e)
336 ENDDO
337C
338 DO n=1,nmax
339#include "vectorize.inc"
340 DO j=1,nindx
341 i=index(j)
342 dpla(i)=dpla_j(i)
343cc PLA_I = PLA(I) + DPLA(I)
344 yld_i = yld(i) + h(i)*dpla(i)
345C
346 dr(i) =half*e*dpla(i)/yld_i
347 pp(i) =one/(one+dr(i)*nu1)
348 qq(i) =one/(one + three*dr(i)*nu2)
349 p2 =pp(i)*pp(i)
350 q2 =qq(i)*qq(i)
351 f =aa(i)*p2+bb(i)*q2-yld_i*yld_i
352 df =-(aa(i)*nu1*p2*pp(i)+three*bb(i)*nu2*q2*qq(i))
353 . *(e-two*dr(i)*h(i))/yld_i
354 . -two*h(i)*yld_i
355 df = sign(max(abs(df),em20),df)
356 IF(dpla(i)>zero) THEN
357 dpla_j(i)=max(zero,dpla(i)-f/df)
358 ELSE
359 dpla_j(i)=zero
360 ENDIF
361 ENDDO
362 ENDDO
363C------------------------------------------
364C CONTRAINTES PLASTIQUEMENT ADMISSIBLES
365C------------------------------------------
366#include "vectorize.inc"
367 DO j=1,nindx
368 i=index(j)
369 pla(i) = pla(i) + dpla(i)
370 s1=(signxx(i)+signyy(i))*pp(i)
371 s2=(signxx(i)-signyy(i))*qq(i)
372 signxx(i)=half*(s1+s2)
373 signyy(i)=half*(s1-s2)
374 signxy(i)=signxy(i)*qq(i)
375 IF (inloc == 0) THEN
376 dezz = - nu3*dr(i)*s1/e
377 thk(i) = thk(i) + dezz*thkly(i)*off(i)
378 ENDIF
379 ENDDO
380 ENDIF
381C------------------------------------------
382 ELSEIF(iflag(1)==2)THEN
383C radial projection on the deviator on a reduced criterion
384C Elastic projection in z => sig33 = 0
385C the coefficient of reduction of the criterion is such that
386C one is on the criterion after the 2 projections
387 DO i=1,nel0
388 p = -(signxx(i)+signyy(i))*third
389 s11 = signxx(i)+p
390 s22 = signyy(i)+p
391C s33 = p = -(S11 + S22)
392 s12 = signxy(i)
393C
394 p2 = p*p
395 vm2= three*(s12*s12 - s11*s22)
396 a = p2*nu4 + vm2
397 vm2= three*p2 + vm2
398 b = p2*nu6
399 c = p2*nu5 - yld(i)*yld(i)
400 r = min(one,(-b + sqrt(max(zero,b*b-a*c)))/max(a ,em20))
401 signxx(i) = s11*r - p
402 signyy(i) = s22*r - p
403 signxy(i) = s12*r
404C signzz = p*r - p
405C proj. signzz = 0.
406 umr = one - r
407 sigz = nnu1*p*umr
408 signxx(i) = signxx(i) + sigz
409 signyy(i) = signyy(i) + sigz
410 svm(i)=sqrt(vm2)
411 dpla(i) = off(i)*svm(i)*umr/g3
412cc DPLA(I) = OFF(I)*SVM(I)*UMR/(G3(I)+H(I))
413 pla(i) = pla(i) + dpla(i)
414 IF (inloc == 0) THEN
415 dezz = dpla(i) * half*(signxx(i)+signyy(i)) / yld(i)
416 dezz=-(depsxx(i)+depsyy(i))*nnu1-nu3*dezz
417 thk(i) = thk(i) + dezz*thkly(i)*off(i)
418 ENDIF
419 IF(r<one) etse(i)= h(i)/(h(i)+e)
420 ENDDO
421 ENDIF
422C calcul de l'effet martensite ----
423 DO i=1,nel0
424 dvm = zero
425 IF(temp(i)<md) THEN
426 dvm = cn*cd*exp((cn-one)*log(max(cd*pla(i),em20)))
427 dvm = dvm*exp(-(cd*pla(i))**cn)
428 dvm = dvm*v0*log(md - temp(i) + one)
429 ENDIF
430
431 vm(i) = vm(i) + max(dvm*dpla(i),zero)
432 vm(i) = min(vm(i), one)
433C the latent heat in the source term.
434 IF(jthe > 0 ) die(i) = (vm(i) - uvar(i,2))*hl
435 uvar(i,2) = vm(i)
436 uvar(i,5) = temp(i)
437 uvar(i,1) = pla(i)
438 ENDDO
439C
440C--------------------------------
441C NON-LOCAL THICKNESS VARIATION
442C--------------------------------
443 IF (inloc > 0) THEN
444 DO i = 1,nel0
445 IF (loff(i) == one) THEN
446 svm(i) = sqrt(signxx(i)*signxx(i) + signyy(i)*signyy(i)
447 . - signxx(i)*signyy(i) + three*signxy(i)*signxy(i))
448 dezz = max(dplanl(i),zero)*half*(signxx(i)+signyy(i))/max(svm(i),em20)
449 dezz = -nu*((signxx(i)-sigoxx(i)+signyy(i)-sigoyy(i))/e) - dezz
450 thk(i) = thk(i) + dezz*thkly(i)*off(i)
451 ENDIF
452 ENDDO
453 ENDIF
454C
455 RETURN
456 END
457C
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
sigeps64c
subroutine sigeps64c(nel0, nuparam, nuvar, mfunc, kfunc, npf, npt0, ipt, iflag, tf, time, timestep, uparam, rho0, area, eint, thkly, epspxx, epspyy, epspxy, epspyz, epspzx, depsxx, depsyy, depsxy, depsyz, depszx, epsxx, epsyy, epsxy, epsyz, epszx, sigoxx, sigoyy, sigoxy, sigoyz, sigozx, signxx, signyy, signxy, signyz, signzx, sigvxx, sigvyy, sigvxy, sigvyz, sigvzx, soundsp, viscmax, thk, pla, uvar, off, ngl, ipm, mat, etse, gs, vol, yld, tempel, die, inloc, dplanl, jthe, loff)
Definition sigeps64c.F:49
vinter
subroutine vinter(tf, iad, ipos, ilen, nel, x, dydx, y)
Definition vinter.F:73