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),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,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,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 . cc,ahs,
163 . bhs, cn, cm, k1,
164 . k2, dh, vm0, temp(mvsiz),
165 . vol0,dvm, vm(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 coef= uparam(24)
202 g3 = three*g
203C
204 nnu1 = nu / (one - nu)
205 nnu2 = nnu1*nnu1
206 nu1 = one/(one-nu)
207 nu2 = one/(one+nu)
208 nu3 = one-nnu1
209 nu4 = one + nnu2 + nnu1
210 nu5 = one + nnu2 - two*nnu1
211 nu6 = half - nnu2 + half*nnu1
212 DO i=1,nel0
213 temp(i) = uparam(22)
214 ENDDO
215C
216 IF (isigi==0) THEN
217 IF(time==0.0)THEN
218 DO i=1,nel0
219 uvar(i,2) = vm0
220 ENDDO
221 ENDIF
222 ENDIF
223C-----------------------------------------------
224 DO i=1,nel0
225C
226 signxx(i)=sigoxx(i) + a1*depsxx(i) + a2*depsyy(i)
227 signyy(i)=sigoyy(i) + a2*depsxx(i) + a1*depsyy(i)
228 signxy(i)=sigoxy(i) + g *depsxy(i)
229 signyz(i)=sigoyz(i) + gs(i)*depsyz(i)
230 signzx(i)=sigozx(i) + gs(i)*depszx(i)
231C
232 soundsp(i) = sqrt(a1/rho0(i))
233 viscmax(i) = zero
234 etse(i) = one
235 ENDDO
236C
237 IF(jthe > 0 ) THEN
238 DO i=1,nel0
239 temp(i) = tempel(i)
240 ENDDO
241 ELSE
242 DO i=1,nel0
243 vol0 = vol(i) * rho0(i)
244 vm(i) = uvar(i,2)
245 temp(i) = temp(i) + (coef*(eint(i,1)+eint(i,2)) + vm(i)*hl) * cp/vol0
246 ENDDO
247 ENDIF
248C
249C compute YLD stress
250 DO i=1,nel0
251 vm(i) = uvar(i,2)
252 aux0 = pla(i)+eps0
253 aux1 = log(max(em20,aux0))
254 aux2 = exp( (cn - one)*aux1 )
255 aux3 = (bhs - ahs) * exp(-cm * aux0 * aux2)
256 kt(i)= k1 + k2*temp(i)
257
258 yld(i) = ( bhs - aux3 ) * kt(i)
259 . + dh*vm(i)
260 h(i) = cm*cn* aux2 * aux3 * kt(i)
261 dpla(i) =zero
262 ENDDO
263C-------------------
264C PROJECTION
265C-------------------
266 IF(iflag(1)==0)THEN
267C radial projection
268 DO i=1,nel0
269 svm(i)=sqrt(signxx(i)*signxx(i)
270 . +signyy(i)*signyy(i)
271 . -signxx(i)*signyy(i)
272 . + three*signxy(i)*signxy(i))
273 r = min(one,yld(i)/max(em20,svm(i)))
274 signxx(i)=signxx(i)*r
275 signyy(i)=signyy(i)*r
276 signxy(i)=signxy(i)*r
277 umr = one - r
278Ctmp DPLA(I) = OFF(I)*SVM(I)*UMR/(G3(I)+H(I))
279 dpla(i) = off(i)*svm(i)*umr/(e)
280 pla(i) = pla(i) + dpla(i)
281 s1=half*(signxx(i)+signyy(i))
282 IF (inloc == 0) THEN
283 dezz = dpla(i) * half*(signxx(i)+signyy(i)) / yld(i)
284 dezz=-(depsxx(i)+depsyy(i))*nnu1-nu3*dezz
285 thk(i) = thk(i) + dezz*thkly(i)*off(i)
286 ENDIF
287 IF(r<1.) etse(i)= h(i)/(h(i)+e)
288 ENDDO
289 ELSEIF(iflag(1)==1)THEN
290C-------------------------
291C von mises criterion
292C-------------------------
293 DO i=1,nel0
294 h(i) = max(zero,h(i))
295 s1=signxx(i)+signyy(i)
296 s2=signxx(i)-signyy(i)
297 s3=signxy(i)
298 aa(i)=fourth*s1*s1
299 bb(i)=three_over_4*s2*s2+three*s3*s3
300 svm(i)=sqrt(aa(i)+bb(i))
301 IF (inloc == 0) THEN
302 dezz = -(depsxx(i)+depsyy(i))*nnu1
303 thk(i) = thk(i) + dezz*thkly(i)*off(i)
304 ENDIF
305 ENDDO
306C-------------------------
307C GATHER PLASTIC FLOW
308C-------------------------
309 nindx=0
310 DO i=1,nel0
311 IF(svm(i)>yld(i).AND.off(i)==1.) THEN
312 nindx=nindx+1
313 index(nindx)=i
314 ENDIF
315 ENDDO
316C
317 IF(nindx/=0) THEN
318C---------------------------
319C DEP EN CONTRAINTE PLANE
320C---------------------------
321 DO j=1,nindx
322 i=index(j)
323 dpla_j(i)=(svm(i)-yld(i))/(g3+h(i))
324 etse(i)= h(i)/(h(i)+e)
325 ENDDO
326C
327 DO n=1,nmax
328#include "vectorize.inc"
329 DO j=1,nindx
330 i=index(j)
331 dpla(i)=dpla_j(i)
332 pla_i = pla(i) + dpla(i)
333C Yld
334 yld_i = ( bhs - (bhs - ahs)*
335 . exp(-cm*exp(cn*log(max(em20,pla_i+eps0))) ))
336 . * kt(i) + dh*vm(i)
337C
338 dr(i) =half*e*dpla(i)/yld_i
339 pp(i) =one/(one+dr(i)*nu1)
340 qq(i) =one/(one+three*dr(i)*nu2)
341 p2 =pp(i)*pp(i)
342 q2 =qq(i)*qq(i)
343 f =aa(i)*p2+bb(i)*q2-yld_i*yld_i
344 df =-(aa(i)*nu1*p2*pp(i)+three*bb(i)*nu2*q2*qq(i))
345 . *(e- two*dr(i)*h(i))/yld_i
346 . - two*h(i)*yld_i
347 df = sign(max(abs(df),em20),df)
348 IF(dpla(i)>zero) THEN
349 dpla_j(i)=max(zero,dpla(i)-f/df)
350 ELSE
351 dpla_j(i)=zero
352 ENDIF
353 ENDDO
354 ENDDO
355C------------------------------------------
356C CONTRAINTES PLASTIQUEMENT ADMISSIBLES
357C------------------------------------------
358#include "vectorize.inc"
359 DO j=1,nindx
360 i=index(j)
361 pla(i) = pla(i) + dpla(i)
362 s1=(signxx(i)+signyy(i))*pp(i)
363 s2=(signxx(i)-signyy(i))*qq(i)
364 signxx(i)=half*(s1+s2)
365 signyy(i)=half*(s1-s2)
366 signxy(i)=signxy(i)*qq(i)
367 IF (inloc == 0) THEN
368 dezz = - nu3*dr(i)*s1/e
369 thk(i) = thk(i) + dezz*thkly(i)*off(i)
370 ENDIF
371 ENDDO
372 ENDIF
373C------------------------------------------
374 ELSEIF(iflag(1)==2)THEN
375C radial projection on the deviator on a reduced criterion
376C Elastic projection in z => sig33 = 0
377C the criterion reduction coefficient is such that
378C one is on the criterion after the 2 projections
379 DO i=1,nel0
380 p = -(signxx(i)+signyy(i))*third
381 s11 = signxx(i)+p
382 s22 = signyy(i)+p
383C s33 = p = -(S11 + S22)
384 s12 = signxy(i)
385C
386 p2 = p*p
387 vm2= three*(s12*s12 - s11*s22)
388 a = p2*nu4 + vm2
389 vm2= three*p2 + vm2
390 b = p2*nu6
391 c = p2*nu5 - yld(i)*yld(i)
392 r = min(one,(-b + sqrt(max(zero,b*b-a*c)))/max(a ,em20))
393 signxx(i) = s11*r - p
394 signyy(i) = s22*r - p
395 signxy(i) = s12*r
396C signzz = p*r - p
397C proj. signzz = 0.
398 umr = one - r
399 sigz = nnu1*p*umr
400 signxx(i) = signxx(i) + sigz
401 signyy(i) = signyy(i) + sigz
402 svm(i)=sqrt(vm2)
403 dpla(i) = off(i)*svm(i)*umr/(three*g)
404ctmp DPLA(I) = OFF(I)*SVM(I)*UMR/(G3(I)+H(I))
405 pla(i) = pla(i) + dpla(i)
406 IF (inloc == 0) THEN
407 dezz = dpla(i) * half*(signxx(i)+signyy(i)) / yld(i)
408 dezz=-(depsxx(i)+depsyy(i))*nnu1-nu3*dezz
409 thk(i) = thk(i) + dezz*thkly(i)*off(i)
410 ENDIF
411 IF(r<one) etse(i)= h(i)/(h(i)+e)
412 ENDDO
413 ENDIF
414C calcul de l'effet martensite ----
415 DO i=1,nel0
416 dvm = half*(one - tanh(cc + cd*temp(i)) )
417
418 argexp= cq/max(em20, temp(i))
419 . + ((cb+one)/cb)
420 . *log( max(em20,one - vm(i))/max(em20,vm(i)))
421 . + pn*log(max(em20,vm(i)))
422
423 dvm = cb*exp( argexp )*dvm/max(ca,em20)
424
425 IF(dpla(i)/=zero) uvar(i,3) = dvm
426 vm(i) = vm(i) + max(dvm*dpla(i),zero)
427 vm(i) = min(vm(i), one)
428C the latent heat in the source term.
429 IF(jthe > 0 ) die(i) = (vm(i) - uvar(i,2))*hl
430 uvar(i,2) = vm(i)
431 uvar(i,4) = temp(i)
432 uvar(i,1) = pla(i)
433 ENDDO
434C--------------------------------
435C NON-LOCAL THICKNESS VARIATION
436C--------------------------------
437 IF (inloc > 0) THEN
438 DO i = 1,nel0
439 IF (loff(i) == one) THEN
440 svm(i) = sqrt(signxx(i)*signxx(i) + signyy(i)*signyy(i)
441 . - signxx(i)*signyy(i) + three*signxy(i)*signxy(i))
442 dezz = max(dplanl(i),zero)*half*(signxx(i)+signyy(i))/max(svm(i),em20)
443 dezz = -nu*((signxx(i)-sigoxx(i)+signyy(i)-sigoyy(i))/e) - dezz
444 thk(i) = thk(i) + dezz*thkly(i)*off(i)
445 ENDIF
446 ENDDO
447 ENDIF
448C
449 RETURN
450 END
451C
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