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