OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps43c.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!|| sigeps43c ../engine/source/materials/mat/mat043/sigeps43c.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 sigeps43c(
32 1 NEL ,NUPARAM ,NUVAR ,NPF ,TF ,
33 2 NFUNC ,IFUNC ,ISRATE ,TIME ,TIMESTEP ,
34 3 UPARAM ,RHO0 ,AREA ,EINT ,THKLY ,
35 4 EPSPXX ,EPSPYY ,EPSPXY ,EPSPYZ ,EPSPZX ,
36 5 DEPSXX ,DEPSYY ,DEPSXY ,DEPSYZ ,DEPSZX ,
37 6 EPSXX ,EPSYY ,EPSXY ,EPSYZ ,EPSZX ,
38 7 SIGOXX ,SIGOYY ,SIGOXY ,SIGOYZ ,SIGOZX ,
39 8 SIGNXX ,SIGNYY ,SIGNXY ,SIGNYZ ,SIGNZX ,
40 9 SOUNDSP ,THK ,PLA ,UVAR ,OFF ,
41 A NGL ,ETSE ,HARDM ,SIGY ,GS ,
42 B DPLA1 ,EPSP ,SHF ,INLOC ,DPLANL ,
43 C SEQ ,LOFF ,ASRATE ,EPSD_PG )
44C=======================================================================
45C Tabulated Hill material model
46C=======================================================================
47C I M P L I C I T T Y P E S
48C-----------------------------------------------
49#include "implicit_f.inc"
50C-----------------------------------------------
51C G L O B A L P A R A M E T E R S
52C-----------------------------------------------
53#include "mvsiz_p.inc"
54C-----------------------------------------------
55C I N P U T A R G U M E N T S
56C-----------------------------------------------
57 INTEGER NEL, NUPARAM, NUVAR, NFUNC ,ISRATE,INLOC
58 INTEGER NGL(NEL),IFUNC(NFUNC)
59 my_real ,INTENT(IN) :: ASRATE
60 my_real ,DIMENSION(NEL) ,INTENT(IN) :: EPSD_PG
61 my_real ,DIMENSION(NEL) ,INTENT(INOUT) :: EPSP
62 my_real
63 . TIME,TIMESTEP,UPARAM(NUPARAM),AREA(NEL),RHO0(NEL),EINT(NEL,2),
64 . THKLY(NEL),PLA(NEL),EPSPXX(NEL),EPSPYY(NEL),GS(NEL),SHF(NEL),
65 . EPSPXY(NEL),EPSPYZ(NEL),EPSPZX(NEL),
66 . DEPSXX(NEL),DEPSYY(NEL),
67 . DEPSXY(NEL),DEPSYZ(NEL),DEPSZX(NEL),
68 . EPSXX(NEL) ,EPSYY(NEL) ,
69 . epsxy(nel) ,epsyz(nel) ,epszx(nel) ,
70 . sigoxx(nel),sigoyy(nel),
71 . sigoxy(nel),sigoyz(nel),sigozx(nel),sigy(nel),hardm(nel),
72 . dplanl(nel),seq(nel)
73 my_real, DIMENSION(NEL), INTENT(IN) :: loff
74C-----------------------------------------------
75C O U T P U T A R G U M E N T S
76C-----------------------------------------------
77 my_real
78 . signxx(nel),signyy(nel),signxy(nel),signyz(nel),signzx(nel),
79 . soundsp(nel),etse(nel)
80C-----------------------------------------------
81C I N P U T O U T P U T A R G U M E N T S
82C-----------------------------------------------
83 my_real uvar(nel,nuvar), off(nel),thk(nel),dpla1(nel)
84C-----------------------------------------------
85C VARIABLES FOR FUNCTION INTERPOLATION
86C-----------------------------------------------
87 INTEGER NPF(*)
88 my_real FINTER ,TF(*)
89C-----------------------------------------------
90C L O C A L V A R I A B L E S
91C-----------------------------------------------
92 INTEGER I,J,NRATE,J1,J2,N,NINDX,NMAX,OPTE,NS,IFUNCE,
93 . IAD1(MVSIZ),IPOS1(MVSIZ),ILEN1(MVSIZ),
94 . IAD2(MVSIZ),IPOS2(MVSIZ),ILEN2(MVSIZ),
95 . jj(mvsiz),index(mvsiz)
96 my_real
97 . a,b,s12,fac,dezz,sigz,s1,s2,s3,dpla,vm2,epst,
98 . a_1,a_2,a_3,p_1,p_2,p_3,pp1,pp2,pp3,dr,dr0,
99 . e(mvsiz),a1(mvsiz),a2(mvsiz),g(mvsiz),g3(mvsiz),
100 . a01,a02,a03,a12,
101 . axx(mvsiz),ayy(mvsiz),axy(mvsiz),a_xy(mvsiz),
102 . dydx1(mvsiz),dydx2(mvsiz),rate(mvsiz,2),svm(mvsiz),
103 . svm0(mvsiz),yld(mvsiz),y1(mvsiz),y2(mvsiz),
104 . yfac(mvsiz,2),nnu1,nu,nu2(mvsiz),
105 . nu3(mvsiz),nu4(mvsiz),nu5(mvsiz),jq(mvsiz),jq2(mvsiz),
106 . dpla_i(mvsiz),dpla_j(mvsiz),fail(mvsiz),epsmax,
107 . epsr1,epsr2,s11(mvsiz),s22(mvsiz),
108 . b_1(mvsiz),b_2(mvsiz),b_3(mvsiz),q12(mvsiz),q21(mvsiz),
109 . err,f,df,pla_i,yld_i,h(mvsiz),
110 . fisokin,hk(mvsiz),fhk,fa01,fa02,fa03,nu1,
111 . ce,einf,escale(mvsiz),c1(mvsiz),normxx,normyy,
112 . dydxe(mvsiz)
113 DATA nmax/2/,ns/10/
114C=======================================================================
115 nrate = nint(uparam(1)) ! nfunc sig(eps)
116 nu = uparam(6)
117 a01 = uparam(7)
118 a02 = uparam(8)
119 a03 = uparam(9)
120 a12 = uparam(10)
121 epsmax=uparam(ns+2*nrate+1)
122 IF(epsmax == zero)THEN
123 IF(tf(npf(ifunc(1)+1)-1) == zero)THEN
124 epsmax=tf(npf(ifunc(1)+1)-2)
125 ELSE
126 epsmax= ep30
127 ENDIF
128 ENDIF
129 nnu1 = nu / (one - nu)
130 epsr1 =uparam(ns+2*nrate+2)
131 epsr2 =uparam(ns+2*nrate+3)
132 fisokin =uparam(ns+2*nrate+8)
133c------------------------------------------------------
134 opte =uparam(ns+2*nrate+10)
135 einf =uparam(ns+2*nrate+11)
136 ce =uparam(ns+2*nrate+12)
137 DO i=1,nel
138 e(i) = uparam(2)
139 a1(i) = uparam(3)
140 a2(i) = uparam(4)
141 g(i) = uparam(5)
142 g3(i) = uparam(ns+2*nrate+5)
143 ENDDO
144 IF (opte == 1)THEN
145 DO i=1,nel
146 IF (pla(i) > zero)THEN
147 ifunce = uparam(ns+2*nrate+9)
148 escale(i) = finter(ifunc(ifunce),pla(i),npf,tf,dydxe(i))
149 e(i) = escale(i)* e(i)
150 a1(i) = e(i)/(one - nu*nu)
151 a2(i) = nu*a1(i)
152 g(i) = half*e(i)/(one+nu)
153 g3(i) = three*g(i)
154 gs(i) = g(i) * shf(i)
155 ENDIF
156 ENDDO
157 ELSEIF ( ce /= zero) THEN
158 DO i=1,nel
159 IF(pla(i) > zero)THEN
160 e(i) = e(i)-(e(i)-einf)*(one-exp(-ce*pla(i)))
161 a1(i) = e(i)/(one - nu*nu)
162 a2(i) = nu*a1(i)
163 g(i) = half*e(i)/(one+nu)
164 g3(i) = three*g(i)
165 gs(i) = g(i) * shf(i)
166 ENDIF
167 ENDDO
168 ENDIF
169c--------------------------------------------------------------
170 dpla1(1:nel) = zero
171 DO i=1,nel
172 signxx(i)=sigoxx(i) - uvar(i,2) +a1(i)*depsxx(i)+a2(i)*depsyy(i)
173 signyy(i)=sigoyy(i) - uvar(i,3) +a2(i)*depsxx(i)+a1(i)*depsyy(i)
174 signxy(i)=sigoxy(i) - uvar(i,4) +g(i) *depsxy(i)
175 signyz(i)=sigoyz(i)+gs(i)*depsyz(i)
176 signzx(i)=sigozx(i)+gs(i)*depszx(i)
177C
178 soundsp(i) = sqrt(a1(i)/rho0(i))
179 etse(i) = one
180C-------------------
181C STRAIN RATE
182C-------------------
183 IF (israte == 0) THEN
184 epsp(i) = half*( abs(epspxx(i)+epspyy(i))
185 . + sqrt( (epspxx(i)-epspyy(i))*(epspxx(i)-epspyy(i))
186 . + epspxy(i)*epspxy(i) ) )
187 ELSE
188 epsp(i) = asrate*epsd_pg(i) + (one-asrate)*epsp(i)
189 ENDIF
190C-------------------
191C STRAIN
192C-------------------
193 epst = half*( epsxx(i)+epsyy(i)
194 . + sqrt( (epsxx(i)-epsyy(i))*(epsxx(i)-epsyy(i))
195 . + epsxy(i)*epsxy(i) ) )
196 fail(i) = max(zero,min(one,(epsr2-epst)/(epsr2-epsr1)))
197C
198 ENDDO
199C-------------------
200C HARDENING LAW
201C-------------------
202 DO i=1,nel
203 jj(i) = 1
204 ENDDO
205 DO i=1,nel
206 DO j=2,nrate-1
207 IF(epsp(i)>=uparam(ns+j)) jj(i) = j
208 ENDDO
209 rate(i,1)=uparam(ns+jj(i))
210 rate(i,2)=uparam(ns+jj(i)+1)
211 yfac(i,1)=uparam(ns+nrate+jj(i))
212 yfac(i,2)=uparam(ns+nrate+jj(i)+1)
213 ENDDO
214 DO i=1,nel
215 j1 = jj(i)
216 j2 = j1+1
217 ipos1(i) = nint(uvar(i,j1+4))
218 iad1(i) = npf(ifunc(j1)) / 2 + 1
219 ilen1(i) = npf(ifunc(j1)+1) / 2 - iad1(i) - ipos1(i)
220 ipos2(i) = nint(uvar(i,j2+4))
221 iad2(i) = npf(ifunc(j2)) / 2 + 1
222 ilen2(i) = npf(ifunc(j2)+1) / 2 - iad2(i) - ipos2(i)
223 ENDDO
224C
225 CALL vinter(tf,iad1,ipos1,ilen1,nel,pla,dydx1,y1)
226 CALL vinter(tf,iad2,ipos2,ilen2,nel,pla,dydx2,y2)
227C
228 DO i=1,nel
229 j1 = jj(i)
230 j2 = j1+1
231 y1(i)=y1(i)*yfac(i,1)
232 y2(i)=y2(i)*yfac(i,2)
233 fac = (epsp(i) - rate(i,1))/(rate(i,2) - rate(i,1))
234 yld(i) = fail(i)*(y1(i) + fac*(y2(i)-y1(i)))
235 h(i) = fail(i)*(dydx1(i) + fac*(dydx2(i)-dydx1(i)))
236 uvar(i,j1+4) = ipos1(i)
237 uvar(i,j2+4) = ipos2(i)
238C ECROUISSAGE CINEMATIQUE
239 y1(i)=tf(npf(ifunc(j1))+1)
240 y2(i)=tf(npf(ifunc(j2))+1)
241 yld(i) = (one-fisokin) * yld(i) +
242 . fisokin * (fail(i)*(y1(i) + fac*(y2(i)-y1(i))))
243 yld(i) = max(yld(i),em20)
244 sigy(i) = yld(i)
245 ENDDO
246C-------------------------
247C HILL CRITERION
248C-------------------------
249 DO i=1,nel
250 h(i) = max(zero,h(i))
251 s1=a01*signxx(i)*signxx(i)
252 s2=a02*signyy(i)*signyy(i)
253 s3=a03*signxx(i)*signyy(i)
254 axy(i)=a12*signxy(i)*signxy(i)
255 svm(i)=sqrt(s1+s2-s3+axy(i))
256 IF (inloc == 0) THEN
257 dezz = -(depsxx(i)+depsyy(i))*nnu1
258 thk(i) = thk(i) + dezz*thkly(i)*off(i)
259 ENDIF
260 ENDDO
261C-------------------------
262C GATHER PLASTIC FLOW
263C-------------------------
264 nindx=0
265 DO i=1,nel
266 IF(svm(i)>yld(i).AND.off(i) == one) THEN
267 nindx=nindx+1
268 index(nindx)=i
269 ENDIF
270 ENDDO
271C---------------------------
272 IF (nindx > 0) THEN
273C---------------------------
274C DEP EN CONTRAINTE PLANE
275C---------------------------
276 DO j=1,nindx
277 i=index(j)
278 dpla_j(i)=(svm(i)-yld(i))/(g3(i)+h(i))
279 etse(i)= h(i)/(h(i)+e(i))
280C +++cette partie peut etre neglige quand FHK est tres petit
281 hk(i) = h(i)*fisokin
282 fhk = four_over_3*hk(i)/a1(i)
283 fa01 =a01*fhk
284 fa02 =a02*fhk
285 fa03 =a03*fhk
286 nu1 =nu+half*fhk
287C ---
288 nu2(i) = 1.-nu1*nu1+fhk*fhk
289 nu3(i) = nu1*half
290 nu4(i) = half*(one-nu)
291 nu5(i) = one-nnu1
292 s1=a01*nu1*two-a03-fa03
293 s2=a02*nu1*two-a03-fa03
294 s12=a03-nu1*(a01+a02)+fa03
295 s3=sqrt(nu2(i)*(a01-a02)*(a01-a02)+s12*s12)
296 IF (abs(s1)<em20) THEN
297 q12(i)=zero
298 ELSE
299 q12(i)=-(a01-a02+s3+fa01-fa02)/s1
300 ENDIF
301 IF (abs(s2)<em20) THEN
302 q21(i)=zero
303 ELSE
304 q21(i)=(a01-a02+s3+fa01-fa02)/s2
305 ENDIF
306 jq(i)=one/(1-q12(i)*q21(i))
307 jq2(i)=jq(i)*jq(i)
308 a=a01*q12(i)
309 b=a02*q21(i)
310 a_1=(a01+a03*q21(i)+b*q21(i))*jq2(i)
311 a_2=(a02+a03*q12(i)+a*q12(i))*jq2(i)
312 a_3=(a+b)*jq2(i)*two+a03*(jq2(i)*two-jq(i))
313 s11(i)=signxx(i)+signyy(i)*q12(i)
314 s22(i)=q21(i)*signxx(i)+signyy(i)
315 axx(i)=a_1*s11(i)*s11(i)
316 ayy(i)=a_2*s22(i)*s22(i)
317 a_xy(i)=a_3*s11(i)*s22(i)
318 a=a03*nu3(i)
319 b=s3*jq(i)
320 b_1(i)=a02-a-b+fa02
321 b_2(i)=a01-a+b+fa01
322 b_3(i)=a12*(nu4(i)+half*fhk)
323 h(i) = h(i)-hk(i)
324 h(i) = max(zero,h(i))
325 ENDDO
326C
327 DO n=1,nmax
328 DO j=1,nindx
329 i=index(j)
330 dpla_i(i)=dpla_j(i)
331 IF(dpla_i(i)>zero) THEN
332 yld_i =yld(i)+h(i)*dpla_i(i)
333 dr =a1(i)*dpla_i(i)/yld_i
334 p_1=one/(one + b_1(i)*dr)
335 pp1=p_1*p_1
336 p_2=one/(one+b_2(i)*dr)
337 pp2=p_2*p_2
338 p_3=one/(one+b_3(i)*dr)
339 pp3=p_3*p_3
340 f =axx(i)*pp1+ayy(i)*pp2-a_xy(i)*p_1*p_2+axy(i)*pp3
341 . -yld_i*yld_i
342 df =-((axx(i)*p_1-a_xy(i)*p_2*half)*pp1*b_1(i)+
343 . (ayy(i)*p_2-a_xy(i)*p_1*half)*pp2*b_2(i)+
344 . axy(i)*pp3*p_3*b_3(i))*(a1(i)-dr*h(i))/yld_i
345 . -h(i)*yld_i
346 dpla_j(i)=max(zero,dpla_i(i)-f*half/df)
347 ELSE
348 dpla_j(i)=zero
349 ENDIF
350 ENDDO
351 ENDDO
352C------------------------------------------
353C CONTRAINTES PLASTIQUEMENT ADMISSIBLES
354C------------------------------------------
355 DO j=1,nindx
356 i=index(j)
357 pla(i) = pla(i) + dpla_j(i)
358 yld_i =yld(i)+h(i)*dpla_j(i)
359 dr0 =dpla_j(i)/yld_i
360 dr =a1(i)*dr0
361 p_1=one/(one + b_1(i)*dr)
362 p_2=one/(one + b_2(i)*dr)
363 p_3=one/(one + b_3(i)*dr)
364 s1=s11(i)*p_1
365 s2=s22(i)*p_2
366 signxx(i)=jq(i)*(s1-s2*q12(i))
367 signyy(i)=jq(i)*(s2-s1*q21(i))
368 signxy(i)=signxy(i)*p_3
369 s1=a01*signxx(i)+a02*signyy(i)
370 . -a03*(signxx(i)+signyy(i))*half
371 IF (inloc == 0) THEN
372 dezz = - nu5(i)*dpla_j(i)*s1/yld_i
373 thk(i) = thk(i) + dezz*thkly(i)*off(i)
374 ENDIF
375 s1 =a03*half
376 p_1=a01*signxx(i)-s1*signyy(i)
377 p_2=a02*signyy(i)-s1*signxx(i)
378 p_3=a12*signxy(i)
379 dr0 = two_third*dr0*hk(i)
380 uvar(i,2) = uvar(i,2) + (two*p_1+p_2)*dr0
381 uvar(i,3) = uvar(i,3) + (two*p_2+p_1)*dr0
382 uvar(i,4) = uvar(i,4) + half*p_3*dr0
383 dpla1(i) = dpla_j(i)
384 ENDDO
385C
386 DO i=1,nel
387 IF(pla(i)>epsmax.AND.off(i) == one)off(i)=four_over_5
388 ENDDO
389 ENDIF
390 DO i=1,nel
391 signxx(i) = signxx(i) + uvar(i,2)
392 signyy(i) = signyy(i) + uvar(i,3)
393 signxy(i) = signxy(i) + uvar(i,4)
394 ENDDO
395C--------------------------------
396C HARDENING MODULUS
397C--------------------------------
398 DO i=1,nel
399 hardm(i) = h(i)
400 ENDDO
401C---------------------------------------------------------------
402C EQUIVALENT STRESS OUTPUT AND NON-LOCAL THICKNESS VARIATION
403C---------------------------------------------------------------
404 DO i = 1,nel
405 s1 = a01*signxx(i)*signxx(i)
406 s2 = a02*signyy(i)*signyy(i)
407 s3 = a03*signxx(i)*signyy(i)
408 axy(i) = a12*signxy(i)*signxy(i)
409 seq(i) = sqrt(s1+s2-s3+axy(i))
410 IF ((inloc > 0).AND.(loff(i) == one)) THEN
411 normxx = (two*a01*signxx(i) - a03*signyy(i))/(max(two*seq(i),em20))
412 normyy = (two*a02*signyy(i) - a03*signxx(i))/(max(two*seq(i),em20))
413 dezz = max(dplanl(i),zero)*(normxx + normyy)
414 dezz = -nu*((signxx(i)-sigoxx(i)+signyy(i)-sigoyy(i))/e(i)) - dezz
415 thk(i) = thk(i) + dezz*thkly(i)*off(i)
416 ENDIF
417 ENDDO
418!-----------
419 RETURN
420 END
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
sigeps43c
subroutine sigeps43c(nel, nuparam, nuvar, npf, tf, nfunc, ifunc, israte, 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, soundsp, thk, pla, uvar, off, ngl, etse, hardm, sigy, gs, dpla1, epsp, shf, inloc, dplanl, seq, loff, asrate, epsd_pg)
Definition sigeps43c.F:44
vinter
subroutine vinter(tf, iad, ipos, ilen, nel, x, dydx, y)
Definition vinter.F:72