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, COEF, 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),COEF(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,NRATE(MVSIZ),J1,J2,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 . epsmax(mvsiz),cc(mvsiz),dd(mvsiz),
168 . cn, cm(mvsiz), k1(mvsiz),
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, eps0(mvsiz),dezz,s1, s2,
176 . s3, dpla_j(mvsiz), yld_i,dr(mvsiz), p2,q2,nnu1,
177 . s11,s22,pla_i,
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 g3 = three*g
205 nnu1 = nu / (one - nu)
206 nnu2 = nnu1*nnu1
207 nu1 = one/(one - nu)
208 nu2 = one/(one + nu)
209 nu3 = one -nnu1
210 nu4 = one + nnu2 + nnu1
211 nu5 = one + nnu2 - two*nnu1
212 nu6 = half - nnu2 + half*nnu1
213 DO i=1,nel0
214 temp(i) = uparam(12)
215C latent heat
216 hl = uparam(15)
217 coef(i) = one
218C
219 ENDDO
220C
221 IF (isigi==0) THEN
222 IF(time==0.0)THEN
223 DO i=1,nel0
224 uvar(i,1)=zero
225 uvar(i,2)=zero
226 uvar(i,3) = zero
227 uvar(i,4) = zero
228 uvar(i,5) = zero
229 uvar(i,6) = zero
230 ENDDO
231 ENDIF
232 ENDIF
233C-----------------------------------------------
234C Calcul de la tempurature. (conduction ou adiabatique)
235C--------------------
236 IF(jthe > 0 ) THEN
237 DO i=1,nel0
238 temp(i) = tempel(i)
239 ENDDO
240 ELSE
241 DO i=1,nel0
242 vol0 = vol(i) * rho0(i)
243 vm(i) = uvar(i,2)
244 temp(i) = temp(i)
245 . + coef(i)*cp*(eint(i,1)+ eint(i,2))/vol0
246 . + cp*vm(i)*hl/vol0
247 ENDDO
248 ENDIF
249C
250 DO i=1,nel0
251C
252 signxx(i)=sigoxx(i) + a1*depsxx(i) + a2*depsyy(i)
253 signyy(i)=sigoyy(i) + a2*depsxx(i) + a1*depsyy(i)
254 signxy(i)=sigoxy(i) + g *depsxy(i)
255 signyz(i)=sigoyz(i) + gs(i)*depsyz(i)
256 signzx(i)=sigozx(i) + gs(i)*depszx(i)
257 vm(i) = uvar(i,2)
258 dpla(i) = zero
259C ...
260
261 soundsp(i) = sqrt(a1/rho0(i))
262 viscmax(i) = zero
263 etse(i) = one
264C
265 ENDDO
266C compute YLD stress
267 DO i=1,nel0
268 ipos1(i) = nint(uvar(i,3))
269 iad1(i) = npf(ifunc(1)) / 2 + 1
270 ilen1(i) = npf(ifunc(1)+1) / 2 - iad1(i) - ipos1(i)
271 ipos2(i) = nint(uvar(i,4))
272 iad2(i) = npf(ifunc(2)) / 2 + 1
273 ilen2(i) = npf(ifunc(2)+1) / 2 - iad2(i) - ipos2(i)
274 ENDDO
275C
276 CALL vinter(tf,iad1,ipos1,ilen1,nel0,pla,dydx1,y1)
277 CALL vinter(tf,iad2,ipos2,ilen2,nel0,pla,dydx2,y2)
278C Yld stress
279 DO i=1,nel0
280 y1(i)=y1(i)*yfac(1)
281 y2(i)=y2(i)*yfac(2)
282 fac = vm(i)/vmc
283 yld(i) = y1(i) + fac*(y2(i)-y1(i))
284 yld(i) = max(yld(i),em20)
285C ecrouissage
286 dydx1(i)=dydx1(i)*yfac(1)
287 dydx2(i)=dydx2(i)*yfac(2)
288 h(i) = (dydx1(i) + fac*(dydx2(i)-dydx1(i)))
289 uvar(i,3) = ipos1(i)
290 uvar(i,4) = ipos2(i)
291 ENDDO
292C-------------------
293C PROJECTION
294C-------------------
295 IF(iflag(1)==0)THEN
296C projection radiale
297 DO i=1,nel0
298 svm(i)=sqrt(signxx(i)*signxx(i)
299 . +signyy(i)*signyy(i)
300 . -signxx(i)*signyy(i)
301 . +three*signxy(i)*signxy(i))
302 r = min(one,yld(i)/max(em20,svm(i)))
303 signxx(i)=signxx(i)*r
304 signyy(i)=signyy(i)*r
305 signxy(i)=signxy(i)*r
306 umr = one - r
307cc DPLA(I) = OFF(I)*SVM(I)*UMR/(G3(I)+H(I))
308 dpla(i) = off(i)*svm(i)*umr/e
309 pla(i) = pla(i) + dpla(i)
310 s1=half*(signxx(i)+signyy(i))
311 IF (inloc == 0) THEN
312 dezz = dpla(i) * half*(signxx(i)+signyy(i)) / yld(i)
313 dezz=-(depsxx(i)+depsyy(i))*nnu1-nu3*dezz
314 thk(i) = thk(i) + dezz*thkly(i)*off(i)
315 ENDIF
316 IF(r<one) etse(i)= h(i)/(h(i)+e)
317 ENDDO
318 ELSEIF(iflag(1)==1)THEN
319C-------------------------
320C CRITERE DE VON MISES
321C-------------------------
322 DO i=1,nel0
323 h(i) = max(zero,h(i))
324 s1=signxx(i)+signyy(i)
325 s2=signxx(i)-signyy(i)
326 s3=signxy(i)
327 aa(i)=fourth*s1*s1
328 bb(i)=three_over_4*s2*s2+three*s3*s3
329 svm(i)=sqrt(aa(i)+bb(i))
330 IF (inloc == 0) THEN
331 dezz = -(depsxx(i)+depsyy(i))*nnu1
332 thk(i) = thk(i) + dezz*thkly(i)*off(i)
333 ENDIF
334 ENDDO
335C-------------------------
336C GATHER PLASTIC FLOW
337C-------------------------
338 nindx=0
339 DO i=1,nel0
340 IF(svm(i)>yld(i).AND.off(i)==one) THEN
341 nindx=nindx+1
342 index(nindx)=i
343 ENDIF
344 ENDDO
345C
346 IF(nindx/=0) THEN
347C---------------------------
348C DEP EN CONTRAINTE PLANE
349C---------------------------
350 DO j=1,nindx
351 i=index(j)
352 dpla_j(i)=(svm(i)-yld(i))/(g3+h(i))
353 etse(i)= h(i)/(h(i)+e)
354 ENDDO
355C
356 DO n=1,nmax
357#include "vectorize.inc"
358 DO j=1,nindx
359 i=index(j)
360 dpla(i)=dpla_j(i)
361cc PLA_I = PLA(I) + DPLA(I)
362 yld_i = yld(i) + h(i)*dpla(i)
363C
364 dr(i) =half*e*dpla(i)/yld_i
365 pp(i) =one/(one+dr(i)*nu1)
366 qq(i) =one/(one + three*dr(i)*nu2)
367 p2 =pp(i)*pp(i)
368 q2 =qq(i)*qq(i)
369 f =aa(i)*p2+bb(i)*q2-yld_i*yld_i
370 df =-(aa(i)*nu1*p2*pp(i)+three*bb(i)*nu2*q2*qq(i))
371 . *(e-two*dr(i)*h(i))/yld_i
372 . -two*h(i)*yld_i
373 df = sign(max(abs(df),em20),df)
374 IF(dpla(i)>zero) THEN
375 dpla_j(i)=max(zero,dpla(i)-f/df)
376 ELSE
377 dpla_j(i)=zero
378 ENDIF
379 ENDDO
380 ENDDO
381C------------------------------------------
382C CONTRAINTES PLASTIQUEMENT ADMISSIBLES
383C------------------------------------------
384#include "vectorize.inc"
385 DO j=1,nindx
386 i=index(j)
387 pla(i) = pla(i) + dpla(i)
388 s1=(signxx(i)+signyy(i))*pp(i)
389 s2=(signxx(i)-signyy(i))*qq(i)
390 signxx(i)=half*(s1+s2)
391 signyy(i)=half*(s1-s2)
392 signxy(i)=signxy(i)*qq(i)
393 IF (inloc == 0) THEN
394 dezz = - nu3*dr(i)*s1/e
395 thk(i) = thk(i) + dezz*thkly(i)*off(i)
396 ENDIF
397 ENDDO
398 ENDIF
399C------------------------------------------
400 ELSEIF(iflag(1)==2)THEN
401C projection radial sur le deviateur sur un critere reduit
402C projection elastique en z => sig33 = 0
403C le coef. de reduction du critere est tel que
404C l'on se trouve sur le critere apres les 2 projections
405 DO i=1,nel0
406 p = -(signxx(i)+signyy(i))*third
407 s11 = signxx(i)+p
408 s22 = signyy(i)+p
409C s33 = p = -(S11 + S22)
410 s12 = signxy(i)
411C
412 p2 = p*p
413 vm2= three*(s12*s12 - s11*s22)
414 a = p2*nu4 + vm2
415 vm2= three*p2 + vm2
416 b = p2*nu6
417 c = p2*nu5 - yld(i)*yld(i)
418 r = min(one,(-b + sqrt(max(zero,b*b-a*c)))/max(a ,em20))
419 signxx(i) = s11*r - p
420 signyy(i) = s22*r - p
421 signxy(i) = s12*r
422C signzz = p*r - p
423C proj. signzz = 0.
424 umr = one - r
425 sigz = nnu1*p*umr
426 signxx(i) = signxx(i) + sigz
427 signyy(i) = signyy(i) + sigz
428 svm(i)=sqrt(vm2)
429 dpla(i) = off(i)*svm(i)*umr/g3
430cc DPLA(I) = OFF(I)*SVM(I)*UMR/(G3(I)+H(I))
431 pla(i) = pla(i) + dpla(i)
432 IF (inloc == 0) THEN
433 dezz = dpla(i) * half*(signxx(i)+signyy(i)) / yld(i)
434 dezz=-(depsxx(i)+depsyy(i))*nnu1-nu3*dezz
435 thk(i) = thk(i) + dezz*thkly(i)*off(i)
436 ENDIF
437 IF(r<one) etse(i)= h(i)/(h(i)+e)
438 ENDDO
439 ENDIF
440C calcul de l'effet martensite ----
441 DO i=1,nel0
442 dvm = zero
443 IF(temp(i)<md) THEN
444 dvm = cn*cd*exp((cn-one)*log(max(cd*pla(i),em20)))
445 dvm = dvm*exp(-(cd*pla(i))**cn)
446 dvm = dvm*v0*log(md - temp(i) + one)
447 ENDIF
448
449 vm(i) = vm(i) + max(dvm*dpla(i),zero)
450 vm(i) = min(vm(i), one)
451C la chaleur latente dans le terme source.
452 IF(jthe > 0 ) die(i) = (vm(i) - uvar(i,2))*hl
453 uvar(i,2) = vm(i)
454 uvar(i,5) = temp(i)
455 uvar(i,1) = pla(i)
456 ENDDO
457C
458C--------------------------------
459C NON-LOCAL THICKNESS VARIATION
460C--------------------------------
461 IF (inloc > 0) THEN
462 DO i = 1,nel0
463 IF (loff(i) == one) THEN
464 svm(i) = sqrt(signxx(i)*signxx(i) + signyy(i)*signyy(i)
465 . - signxx(i)*signyy(i) + three*signxy(i)*signxy(i))
466 dezz = max(dplanl(i),zero)*half*(signxx(i)+signyy(i))/max(svm(i),em20)
467 dezz = -nu*((signxx(i)-sigoxx(i)+signyy(i)-sigoyy(i))/e) - dezz
468 thk(i) = thk(i) + dezz*thkly(i)*off(i)
469 ENDIF
470 ENDDO
471 ENDIF
472C
473 RETURN
474 END
475C
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, coef, inloc, dplanl, jthe, loff)
Definition sigeps64c.F:49
vinter
subroutine vinter(tf, iad, ipos, ilen, nel, x, dydx, y)
Definition vinter.F:72