OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps60g.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!|| sigeps60g ../engine/source/materials/mat/mat060/sigeps60g.F
25!||--- called by ------------------------------------------------------
26!|| mulawglc ../engine/source/materials/mat_share/mulawglc.F
27!||--- calls -----------------------------------------------------
28!|| finter ../engine/source/tools/curve/finter.F
29!|| inter_rat ../engine/source/materials/mat/mat060/sigeps60c.F
30!|| vinter ../engine/source/tools/curve/vinter.F
31!||====================================================================
32 SUBROUTINE sigeps60g(
33 1 NEL ,NUVAR ,NPF ,
34 2 TF ,TIME ,UPARAM ,RHO0 ,
35 3 AREA ,EINT ,THK0 ,
36 4 EPSPXX ,EPSPYY ,EPSPXY ,EPSPYZ ,EPSPZX ,
37 5 DEPSXX ,DEPSYY ,DEPSXY ,DEPSYZ ,DEPSZX ,
38 5 DEPBXX ,DEPBYY ,DEPBXY ,
39 6 EPSXX ,EPSYY ,EPSXY ,EPSYZ ,EPSZX ,
40 7 SIGOXX ,SIGOYY ,SIGOXY ,SIGOYZ ,SIGOZX ,
41 7 MOMOXX ,MOMOYY ,MOMOXY ,
42 8 SIGNXX ,SIGNYY ,SIGNXY ,SIGNYZ ,SIGNZX ,
43 8 MOMNXX ,MOMNYY ,MOMNXY ,
44 9 SIGVXX ,SIGVYY ,SIGVXY ,SIGVYZ ,SIGVZX ,
45 A SOUNDSP,VISCMAX,THK ,PLA ,UVAR ,
46 B OFF ,NGL ,IPM ,MAT ,ETSE ,
47 C GS ,YLD ,EPSP ,ISRATE ,IPLA ,
48 D SHF )
49C-----------------------------------------------------------------------------
50C ZENG 15/07/97
51C N,M SONT REMPLACES PAR N/H,M/H^2
52C LES DIFFERENTES OPTIONS SONT SUPRIMEES
53C IL NE RESTE QUE : COUPLAGE M-N AVEC R (GAMA DANS LE PROGRAMME)VARIABLE
54C MODIF POUVANT ETRE ENVISAGE: M(X) DONNE PAR USER -> CALCULER R
55C-----------------------------------------------------------------------------
56C I M P L I C I T T Y P E S
57C-----------------------------------------------
58#include "implicit_f.inc"
59C-----------------------------------------------
60C G L O B A L P A R A M E T E R S
61C-----------------------------------------------
62#include "mvsiz_p.inc"
63#include "param_c.inc"
64#include "parit_c.inc"
65#include "com01_c.inc"
66#include "scr03_c.inc"
67C-----------------------------------------------
68C I N P U T A R G U M E N T S
69C-----------------------------------------------
70 INTEGER NEL, NUVAR, NGL(NEL), MAT(NEL),ISRATE,IPLA
71 . ,IPM(NPROPMI,*)
72 my_real
73 . TIME,UPARAM(*),
74 . AREA(NEL),RHO0(NEL),EINT(NEL,2),
75 . THK0(NEL),PLA(NEL),
76 . EPSPXX(NEL),EPSPYY(NEL),
77 . EPSPXY(NEL),EPSPYZ(NEL),EPSPZX(NEL),
78 . DEPSXX(NEL),DEPSYY(NEL),DEPSXY(NEL),
79 . DEPBXX(NEL),DEPBYY(NEL),DEPBXY(NEL),
80 . DEPSYZ(NEL),DEPSZX(NEL),
81 . EPSXX(NEL) ,EPSYY(NEL) ,
82 . EPSXY(NEL) ,EPSYZ(NEL) ,EPSZX(NEL) ,
83 . SIGOXX(NEL),SIGOYY(NEL),SIGOXY(NEL),
84 . MOMOXX(NEL),MOMOYY(NEL),MOMOXY(NEL),
85 . SIGOYZ(NEL),SIGOZX(NEL),
86 . gs(*),epsp(nel),shf(nel)
87C-----------------------------------------------
88C O U T P U T A R G U M E N T S
89C-----------------------------------------------
90 my_real
91 . signxx(nel),signyy(nel),signxy(nel),
92 . momnxx(nel),momnyy(nel),momnxy(nel),
93 . signyz(nel),signzx(nel),
94 . sigvxx(nel),sigvyy(nel),
95 . sigvxy(nel),sigvyz(nel),sigvzx(nel),
96 . soundsp(nel),viscmax(nel),etse(nel)
97C-----------------------------------------------
98C I N P U T O U T P U T A R G U M E N T S
99C-----------------------------------------------
100 my_real
101 . uvar(nel,nuvar), off(nel),thk(nel),yld(nel)
102C-----------------------------------------------
103C VARIABLES FOR FUNCTION INTERPOLATION
104C-----------------------------------------------
105 INTEGER NPF(*)
106 my_real
107 . tf(*),finter
108C Y = FINTER(IFUNC(J),X,NPF,TF,DYDX)
109C Y : Y = F(X)
110C X : X
111C DYDX : F'(X) = DY/DX
112C IFUNC(J): FUNCTION INDEX
113C J : FIRST(J=1), SECOND(J=2) .. FUNCTION USED FOR THIS LAW
114C NPF,TF : FUNCTION PARAMETER
115C-----------------------------------------------
116C L O C A L V A R I A B L E S
117C-----------------------------------------------
118 INTEGER I,J,NRATE,N,NINDX,NMAX,NFUNC,
119 . IAD1(MVSIZ),IPOS1(MVSIZ),ILEN1(MVSIZ),
120 . IAD2(MVSIZ),IPOS2(MVSIZ),ILEN2(MVSIZ),IFUNCE,
121 . jj(mvsiz),index(mvsiz),ifunc(100),nfuncv,
122 . nfuncm, nratem, iadbufv, nn, ipos3(mvsiz),
123 . ipos4(mvsiz),iad3(mvsiz),ilen3(mvsiz),iad4(mvsiz),
124 . ilen4(mvsiz),j1, j2, j3, j4,opte,mx
125 my_real
126 . e(mvsiz),nu,a,b,fac,dezz,s1(mvsiz),s2(mvsiz),
127 . dpla,epst,a1(mvsiz),a2(mvsiz),g(mvsiz),g3(mvsiz),
128 . dydx1(mvsiz),dydx2(mvsiz),rate(mvsiz,4),svm(mvsiz),
129 . y1(mvsiz),y2(mvsiz),dr,
130 . yfac(mvsiz,4),nnu1,nu1(mvsiz),
131 . nu2(mvsiz),nu3(mvsiz),dpla_i,dpla_j(mvsiz),h(mvsiz),
132 . fail(mvsiz),epsmax,epsr1,epsr2,
133 . err,f,df,yld_i,
134 . c1,cp1,cq1,cp2,cq2,sm1(mvsiz),sm2(mvsiz),sm3,fn,fm,fnm,
135 . pn2,qn2,pm2,qm2,dfn,dfm,dfnm,da,db,a_i,b_i,
136 . dfnp,dfnq,dfmp,dfmq,dfnmp,dfnmq,xp,xq,xpg,xqg,
137 . gm(mvsiz),cm(mvsiz),p_m,qm,pnm1,pnm2,qtier(mvsiz),
138 . cnm(mvsiz),am(mvsiz),bm(mvsiz),anm(mvsiz),bnm(mvsiz),
139 . num1(mvsiz),num2(mvsiz),an(mvsiz),bn(mvsiz),
140 . nvm(mvsiz),mvm(mvsiz),nmvm(mvsiz),pn,qn,sn1,sn2,s,
141 . qnm1,qnm2,fnp,fnq,fmp,fmq,fnmp,fnmq,s3,aa,bb,m1,m2,
142 . lfn(mvsiz),qfn(mvsiz),qfnm(mvsiz),rr(mvsiz),
143 . d1,d2,dwt,dwe,dwp,aaa,bbb,ccc,fs,ms,
144 . am1(mvsiz),am2(mvsiz),gama(mvsiz),gama2(mvsiz),
145 . y11, y22, y33, y44,x,y,yp,escale(mvsiz),
146 . x1, x2, x3, x4,y3(mvsiz),y4(mvsiz),dydxe(mvsiz),
147 . dydx3(mvsiz),dydx4(mvsiz),einf,ce
148
149C
150 DATA nmax/2/
151C-----------------------------------------------
152C USER VARIABLES INITIALIZATION
153C-----------------------------------------------
154 IF (ivector/=1) THEN
155 mx = mat(1)
156 nfunc = ipm(10,mx)
157 DO j=1,nfunc
158 ifunc(j)=ipm(10+j,mx)
159 ENDDO
160 iadbufv = ipm(7,mx) - 1
161 nu = uparam(iadbufv+6)
162 nrate = nint(uparam(iadbufv+1))
163 epsmax=uparam(iadbufv+6+2*nrate+1)
164C
165 IF(epsmax==0.)THEN
166 IF(tf(npf(ifunc(1)+1)-1)==zero)THEN
167 epsmax=tf(npf(ifunc(1)+1)-2)
168 ELSE
169 epsmax= ep30
170 ENDIF
171 ENDIF
172C
173 nnu1 = nu / (one - nu)
174 epsr1 =uparam(iadbufv+6+2*nrate+2)
175 IF(epsr1==0.)epsr1=ep30
176 epsr2 =uparam(iadbufv+6+2*nrate+3)
177 IF(epsr2==0.)epsr2=twoep30
178
179c ---
180 opte = uparam(iadbufv+2*nrate+18)
181 einf = uparam(iadbufv+2*nrate+19)
182 ce = uparam(iadbufv+2*nrate+20)
183 DO i=1,nel
184 e(i) = uparam(iadbufv+2)
185 a1(i) = uparam(iadbufv+3)
186 a2(i) = uparam(iadbufv+4)
187 g(i) = uparam(iadbufv+5)
188 g3(i) = 3.*g(i)
189C GCT(I) = G(I)*R
190 c1=thk0(i)*one_over_12
191 am1(i) = a1(i)*c1
192 am2(i) = a2(i)*c1
193 gm(i) = g(i)*c1
194c ---
195 ENDDO
196 IF (opte == 1)THEN
197 DO i=1,nel
198 IF(pla(i) > zero)THEN
199 ifunce = uparam(iadbufv+2*nrate+17)
200 escale(i) =finter(ifunc(ifunce),pla(i),npf,tf,dydxe(i))
201 e(i) = escale(i)* e(i)
202 g(i) = half*e(i)/(one+nu)
203 gs(i) = g(i)*shf(i)
204 g3(i) = three*g(i)
205 a1(i) = e(i)/(one - nu*nu)
206 a2(i) = nu*a1(i)
207 am1(i) = a1(i)*c1
208 am2(i) = a2(i)*c1
209 gm(i) = g(i)*c1
210 ENDIF
211 ENDDO
212 ELSEIF ( ce /= zero) THEN
213 DO i=1,nel
214 IF(pla(i) > zero)THEN
215 e(i) = e(i)-(e(i)-einf)*(one-exp(-ce*pla(i)))
216 g(i) = half*e(i)/(one+nu)
217 gs(i) = g(i)*shf(i)
218 g3(i) = three*g(i)
219 a1(i) = e(i)/(one - nu*nu)
220 a2(i) = nu*a1(i)
221 am1(i) = a1(i)*c1
222 am2(i) = a2(i)*c1
223 gm(i) = g(i)*c1
224 ENDIF
225 ENDDO
226 ENDIF
227 ELSE ! ivector = 1 only
228 nfuncm = 0
229 mx = mat(1)
230 nfuncv = ipm(10,mx)
231 nfuncm = max(nfuncm,nfuncv)
232 DO j=1,nfuncm
233 IF(nfuncv>=j) THEN
234 ifunc(j)=ipm(10+j,mx)
235 ENDIF
236 ENDDO
237 nratem = 0
238C
239 iadbufv = ipm(7,mx) - 1
240 nu = uparam(iadbufv+6)
241 nrate = nint(uparam(iadbufv+1))
242 nratem = max(nratem,nrate)
243 epsmax=uparam(iadbufv+6+2*nrate+1)
244C
245 IF(epsmax==zero)THEN
246 IF(tf(npf(ifunc(1)+1)-1)==zero)THEN
247 epsmax=tf(npf(ifunc(1)+1)-2)
248 ELSE
249 epsmax= 1.e30
250 ENDIF
251 ENDIF
252C
253 nnu1 = nu / (1. - nu)
254 epsr1 =uparam(iadbufv+6+2*nrate+2)
255 IF(epsr1==zero)epsr1=ep30
256 epsr2 =uparam(iadbufv+6+2*nrate+3)
257 IF(epsr2==zero)epsr2=ep30
258c ---
259 opte = uparam(iadbufv+2*nrate+18)
260 einf = uparam(iadbufv+2*nrate+19)
261 ce = uparam(iadbufv+2*nrate+20)
262 DO i=1,nel
263C
264 e(i) = uparam(iadbufv+2)
265 a1(i) = uparam(iadbufv+3)
266 a2(i) = uparam(iadbufv+4)
267 g(i) = uparam(iadbufv+5)
268
269 g3(i) = 3.*g(i)
270C GCT(I) = G(I)*R
271 c1=thk0(i)*one_over_12
272 am1(i) = a1(i)*c1
273 am2(i) = a2(i)*c1
274 gm(i) = g(i)*c1
275c ---
276 ENDDO
277 IF (opte == 1)THEN
278 DO i=1,nel
279 IF(pla(i) > zero)THEN
280 ifunce = uparam(iadbufv+2*nrate+17)
281 escale(i)=finter(ifunc(ifunce),pla(i),npf,tf,dydxe(i))
282 e(i) = escale(i)* e(i)
283 g(i) = half*e(i)/(one+nu)
284 gs(i) = g(i)*shf(i)
285 g3(i) = three*g(i)
286 a1(i) = e(i)/(one - nu*nu)
287 a2(i) = nu*a1(i)
288 am1(i) = a1(i)*c1
289 am2(i) = a2(i)*c1
290 gm(i) = g(i)*c1
291 ENDIF
292 ENDDO
293 ELSEIF ( ce /= zero) THEN
294 DO i=1,nel
295 IF(pla(i) > zero)THEN
296 e(i) = e(i)-(e(i)-einf)*(one-exp(-ce*pla(i)))
297 g(i) = half*e(i)/(one+nu)
298 gs(i) = g(i)*shf(i)
299 g3(i) = three*g(i)
300 a1(i) = e(i)/(one - nu*nu)
301 a2(i) = nu*a1(i)
302 am1(i) = a1(i)*c1
303 am2(i) = a2(i)*c1
304 gm(i) = g(i)*c1
305 ENDIF
306 ENDDO
307 ENDIF
308 ENDIF
309C
310C
311C
312 IF (isigi==0) THEN
313 IF(time==0.0)THEN
314 DO i=1,nel
315 uvar(i,1)=zero
316 uvar(i,2)=zero
317 DO j=1,nrate
318 uvar(i,j+2)=0
319 ENDDO
320 ENDDO
321 ENDIF
322 ENDIF
323C-----------------------------------------------
324C
325 DO i=1,nel
326 signxx(i)=sigoxx(i)+a1(i)*depsxx(i)+a2(i)*depsyy(i)
327 signyy(i)=sigoyy(i)+a2(i)*depsxx(i)+a1(i)*depsyy(i)
328 signxy(i)=sigoxy(i)+g(i) *depsxy(i)
329 momnxx(i)=momoxx(i)+am1(i)*depbxx(i)+am2(i)*depbyy(i)
330 momnyy(i)=momoyy(i)+am2(i)*depbxx(i)+am1(i)*depbyy(i)
331 momnxy(i)=momoxy(i)+gm(i) *depbxy(i)
332 signyz(i)=sigoyz(i)+gs(i) *depsyz(i)
333 signzx(i)=sigozx(i)+gs(i) *depszx(i)
334C
335 soundsp(i) = sqrt(a1(i)/rho0(i))
336 etse(i) = one
337C-------------------
338C STRAIN RATE
339C-------------------
340 IF (israte == 0) epsp(i) = half*( abs(epspxx(i)+epspyy(i))
341 . + sqrt( (epspxx(i)-epspyy(i))*(epspxx(i)-epspyy(i))
342 . + epspxy(i)*epspxy(i) ) )
343C-------------------
344C STRAIN
345C-------------------
346 epst = half*( epsxx(i)+epsyy(i)
347 . + sqrt( (epsxx(i)-epsyy(i))*(epsxx(i)-epsyy(i))
348 . + epsxy(i)*epsxy(i) ) )
349 fail(i) = max(zero,min(one,(epsr2-epst)/(epsr2-epsr1)))
350 ENDDO
351C-------------------
352C CRITERE
353C-------------------
354 DO i=1,nel
355 jj(i) = 1
356 ENDDO
357 IF (ivector/=1) THEN
358C
359 DO i=1,nel
360 DO j=2,nrate-1
361 IF(epsp(i)>=uparam(iadbufv+6+j)) jj(i) = j
362 ENDDO
363 ENDDO
364C
365 ELSE
366 DO j=2,nratem-1
367 DO i=1,nel
368 IF(nrate-1>=j) THEN
369 IF(epsp(i)>=uparam(iadbufv+6+j)) jj(i) = j
370 ENDIF
371 ENDDO
372 ENDDO
373 ENDIF
374C
375 DO i=1,nel
376 IF(jj(i)==1)THEN
377 j1 = jj(i)
378 j2 = j1+1
379 j3 = j2+1
380 j4 = j3+1
381 ELSEIF(jj(i)==(nrate-1))THEN
382 j1 = jj(i) - 2
383 j2 = j1+1
384 j3 = j2+1
385 j4 = j3+1
386 ELSE
387 j1 = jj(i) - 1
388 j2 = j1+1
389 j3 = j2+1
390 j4 = j3+1
391 ENDIF
392 rate(i,1)=uparam(iadbufv + 6 + j1 )
393 rate(i,2)=uparam(iadbufv + 6 + j2 )
394 rate(i,3)=uparam(iadbufv + 6 + j3 )
395 rate(i,4)=uparam(iadbufv + 6 + j4 )
396 yfac(i,1)=uparam(iadbufv+6+nrate + j1 )
397 yfac(i,2)=uparam(iadbufv+6+nrate + j2 )
398 yfac(i,3)=uparam(iadbufv+6+nrate + j3 )
399 yfac(i,4)=uparam(iadbufv+6+nrate + j4 )
400 ENDDO
401C
402C
403 DO i=1,nel
404 IF(jj(i)==1)THEN
405 j1 = jj(i)
406 j2 = j1+1
407 j3 = j2+1
408 j4 = j3+1
409 ELSEIF(jj(i)==(nrate-1))THEN
410 j1 = jj(i) - 2
411 j2 = j1+1
412 j3 = j2+1
413 j4 = j3+1
414 ELSE
415 j1 = jj(i) - 1
416 j2 = j1+1
417 j3 = j2+1
418 j4 = j3+1
419 ENDIF
420 ipos1(i) = nint(uvar(i,j1))
421 iad1(i) = npf(ifunc(j1)) / 2 + 1
422 ilen1(i) = npf(ifunc(j1)+1) / 2 - iad1(i) - ipos1(i)
423 ipos2(i) = nint(uvar(i,j2))
424 iad2(i) = npf(ifunc(j2)) / 2 + 1
425 ilen2(i) = npf(ifunc(j2)+1) / 2 - iad2(i) - ipos2(i)
426 ipos3(i) = nint(uvar(i,j3+4))
427 iad3(i) = npf(ifunc(j3)) / 2 + 1
428 ilen3(i) = npf(ifunc(j3)+1) / 2 - iad3(i) - ipos3(i)
429 ipos4(i) = nint(uvar(i,j4+4))
430 iad4(i) = npf(ifunc(j4)) / 2 + 1
431 ilen4(i) = npf(ifunc(j4)+1) / 2 - iad4(i) - ipos4(i)
432 ENDDO
433C
434 CALL vinter(tf,iad1,ipos1,ilen1,nel,pla,dydx1,y1)
435 CALL vinter(tf,iad2,ipos2,ilen2,nel,pla,dydx2,y2)
436 CALL vinter(tf,iad3,ipos3,ilen3,nel,pla,dydx3,y3)
437 CALL vinter(tf,iad4,ipos4,ilen4,nel,pla,dydx4,y4)
438C
439 DO i=1,nel
440 IF(jj(i)==1)THEN
441 j1 = jj(i)
442 j2 = j1+1
443 j3 = j2+1
444 j4 = j3+1
445 dydx1(i)= dydx1(i)*yfac(i,1)
446 dydx2 (i) = dydx2(i)*yfac(i,2)
447 fac = (epsp(i) - rate(i,1))/(rate(i,2) - rate(i,1))
448 h(i) = fail(i)*(dydx1(i) + fac*(dydx2(i)-dydx1(i)))
449 ELSEIF(jj(i)==(nrate-1))THEN
450 j1 = jj(i) - 2
451 j2 = j1+1
452 j3 = j2+1
453 j4 = j3+1
454 dydx3(i) = dydx3(i)*yfac(i,3)
455 dydx4(i) = dydx4(i)*yfac(i,4)
456 fac = (epsp(i) - rate(i,3))/(rate(i,4) - rate(i,3))
457 h(i) = fail(i)*(dydx3(i) + fac*(dydx4(i)-dydx3(i)))
458 ELSE
459 j1 = jj(i) - 1
460 j2 = j1+1
461 j3 = j2+1
462 j4 = j3+1
463 dydx2(i) = dydx2(i)*yfac(i,2)
464 dydx3(i) = dydx3(i)*yfac(i,3)
465 fac = (epsp(i) - rate(i,2))/(rate(i,3) - rate(i,2))
466 h(i) = fail(i)*(dydx2(i) + fac*(dydx3(i)-dydx2(i)))
467 ENDIF
468 nn = nrate
469 x1 = rate(i,1)
470 x2 = rate(i,2)
471 x3 = rate(i,3)
472 x4 = rate(i,4)
473 x = epsp(i)
474 y11 = y1(i)*yfac(i,1)
475 y22 = y2(i)*yfac(i,2)
476 y33 = y3(i)*yfac(i,3)
477 y44 = y4(i)*yfac(i,4)
478 CALL inter_rat(x1,x2,x3,x4,y11,y22,y33,y44,
479 . x,y,yp,jj(i),nn)
480 yld(i) = fail(i)*y
481 yld(i) = max(yld(i),em20)
482Cc Y11 = DYDX1(I)*YFAC(I,1)
483Cc Y22 = DYDX2(I)*YFAC(I,2)
484Cc Y33 = DYDX3(I)*YFAC(I,3)
485Cc Y44 = DYDX4(I)*YFAC(I,4)
486Cc CALL INTER_RAT(X1,X2,X3,X4,Y11,Y22,Y33,Y44,
487Cc . X,Y,YP,JJ(I),NN)
488Cc H(I) = FAIL(I)*Y
489 uvar(i,j1) = ipos1(i)
490 uvar(i,j2) = ipos2(i)
491 uvar(i,j3) = ipos3(i)
492 uvar(i,j4) = ipos4(i)
493 ENDDO
494 IF(ipla==0) THEN
495C---------------------------
496C RADIAL RETURN
497C---------------------------
498 DO i=1,nel
499 ms=momnxx(i)+momnyy(i)
500 fs=signxx(i)+signyy(i)
501 svm(i) = sixteen*(ms*ms +three*(momnxy(i)*momnxy(i)
502 . - momnxx(i)*momnyy(i)))
503 . + fs*fs+ three*(signxy(i)*signxy(i)-signxx(i)*signyy(i))
504 svm(i) = sqrt(max(svm(i),em20))
505 rr(i) = min(one,yld(i)/svm(i))
506 IF(rr(i)<one) etse(i)= h(i)/(h(i)+e(i))
507 ENDDO
508 DO i=1,nel
509 signxx(i) = signxx(i)*rr(i)
510 signyy(i) = signyy(i)*rr(i)
511 signxy(i) = signxy(i)*rr(i)
512 momnxx(i) = momnxx(i)*rr(i)
513 momnyy(i) = momnyy(i)*rr(i)
514 momnxy(i) = momnxy(i)*rr(i)
515 d1 = signxx(i)-sigoxx(i)
516 d2 = signyy(i)-sigoyy(i)
517 nu = uparam(iadbufv+6)
518 dwe =((signxx(i)+sigoxx(i))*(d1-nu*d2)+
519 . (signyy(i)+sigoyy(i))*(-nu*d1+d2))/e(i)+
520 . (signxy(i)+sigoxy(i))*(signxy(i)-sigoxy(i))/g(i)
521 d1 = momnxx(i)-momoxx(i)
522 d2 = momnyy(i)-momoyy(i)
523 dwe =dwe+ twelve*(
524 . ((momnxx(i)+momoxx(i))*(d1-nu*d2)
525 . +(momnyy(i)+momoyy(i))*(-nu*d1+d2))/e(i)
526 . +(momnxy(i)+momoxy(i))*(momnxy(i)-momoxy(i))/g(i) )
527 dwt = (signxx(i)+sigoxx(i))*depsxx(i)+
528 . (signyy(i)+sigoyy(i))*depsyy(i)+
529 . (signxy(i)+sigoxy(i))*depsxy(i)
530 dwt = dwt+thk0(i)*((momnxx(i)+momoxx(i))*depbxx(i)+
531 . (momnyy(i)+momoyy(i))*depbyy(i)+
532 . (momnxy(i)+momoxy(i))*depbxy(i))
533 dwp =dwt-dwe
534 dpla = off(i)* max(zero,half*dwp/yld(i))
535 pla(i)=pla(i) + dpla
536 aaa = abs(dwe)
537 bbb = max(zero,dwp)
538 ccc = max(em20,aaa+bbb)
539 dezz = - (depsxx(i)+depsyy(i)) * (nnu1*aaa + bbb) / ccc
540 thk(i) = thk(i) * (one + dezz*off(i))
541 ENDDO
542C
543 ELSEIF(codvers<44)THEN
544C IF(CODVERS<44.AND.IPLA==1) THEN
545C-------------------------
546C CRITERE DE VON MISES
547C-------------------------
548 DO i=1,nel
549C-------------------------------------------------------------------------
550C GAMA (L'INVERSE DE GAMA DANS LA FORMULE) A ETE PRIS A 2/3 PAR DEFAUT
551C-------------------------------------------------------------------------
552 c1 = pla(i)*e(i)
553 gama(i) = three_half*(c1+yld(i))/(three_half*c1+yld(i))
554 gama2(i) = gama(i)*gama(i)
555 cm(i) = sixteen*gama2(i)
556 cnm(i) = sqr16_3*gama(i)
557 qtier(i) = four_over_3*gama2(i)
558 h(i) = max(zero,h(i))
559 s1(i) = (signxx(i)+signyy(i))*half
560 s2(i) = (signxx(i)-signyy(i))*half
561 s3 = signxy(i)
562 sm1(i) = (momnxx(i)+momnyy(i))*half
563 sm2(i) = (momnxx(i)-momnyy(i))*half
564 sm3 = momnxy(i)
565 an(i) = s1(i)*s1(i)
566 bn(i) = three*(s2(i)*s2(i)+s3*s3)
567 nvm(i) = an(i)+bn(i)
568 am(i) = sm1(i)*sm1(i)*cm(i)
569 bm(i) =three*(sm2(i)*sm2(i)+sm3*sm3)*cm(i)
570 mvm(i) = am(i)+bm(i)
571 anm(i) = s1(i)*sm1(i)*cnm(i)
572 bnm(i) = three*(s2(i)*sm2(i)+s3*sm3)*cnm(i)
573 nmvm(i) = anm(i)+bnm(i)
574 svm(i) = sqrt(nvm(i)+mvm(i)+abs(nmvm(i)))
575 dezz = -(depsxx(i)+depsyy(i))*nnu1
576 thk(i) = thk(i) +thk(i)* dezz*off(i)
577 ENDDO
578C-------------------------
579C GATHER PLASTIC FLOW
580C-------------------------
581 nindx=0
582 DO i=1,nel
583 IF(svm(i)>yld(i).AND.off(i)==one) THEN
584 nindx=nindx+1
585 index(nindx)=i
586 ENDIF
587 ENDDO
588 IF(nindx==0) RETURN
589C---------------------------
590C DEP EN CONTRAINTE PLANE
591C---------------------------
592 DO j=1,nindx
593 i=index(j)
594 nu=uparam(iadbufv+6)
595 nu1(i) = one/(one - nu)
596 nu2(i) = one/(one + nu)
597 nu3(i) = one - nnu1
598 num1(i) = nu1(i)*qtier(i)
599 num2(i) = nu2(i)*qtier(i)
600 dpla_j(i)=(svm(i)-yld(i))/(g3(i)*qtier(i)+h(i))
601 etse(i)= h(i)/(h(i)+e(i))
602 ENDDO
603C-------------------------------
604C TIENT COMPTE DU COUPLAGE
605C-------------------------------
606 DO n=1,nmax
607#include "vectorize.inc"
608 DO j=1,nindx
609 i=index(j)
610 dpla_i=dpla_j(i)
611 yld_i =yld(i)+h(i)*dpla_i
612 dr =half*e(i)*dpla_i/yld_i
613 xp =dr*nu1(i)
614 xq =three*dr*nu2(i)
615 xpg =xp*zep444*gama2(i)
616 xqg =xq*zep444*gama2(i)
617 c1=one+qtier(i)
618 da=c1+twop444*gama2(i)*xp
619 db=c1+twop444*gama2(i)*xq
620 a=one +(da+c1)*xp*half
621 b=one +(db+c1)*xq*half
622 a_i=one/a
623 b_i=one/b
624 aa=a_i*a_i
625 bb=b_i*b_i
626 dfnp=fivep5+sixteenp5*xpg
627 fnp=one+(dfnp+fivep5)*xpg*half
628 dfnq=fivep5+sixteenp5*xqg
629 fnq=one+(dfnq+fivep5)*xqg*half
630 dfmp=onep8333*(xp+one)
631 fmp=one+(dfmp+onep8333)*xp*half
632 dfmq=onep8333*(xq+one)
633 fmq=one+(dfmq+onep8333)*xq*half
634 dfnmp=-twop444*xp*gama2(i)
635 fnmp=one+dfnmp*xp*half
636 dfnmq=-twop444*xq*gama2(i)
637 fnmq=one+dfnmq*xq*half
638 fn=aa*fnp*an(i)+bb*fnq*bn(i)
639 fm=aa*fmp*am(i)+bb*fmq*bm(i)
640 fnm=aa*fnmp*anm(i)+bb*fnmq*bnm(i)
641 IF (fnm<zero) THEN
642 fnm=-fnm
643 s=-one
644 ELSE
645 s=one
646 ENDIF
647 f =fn+fm+fnm-yld_i*yld_i
648 c1=nu1(i)*aa*a_i
649 cp1=c1*a
650 cp2=c1*da*two
651 c1=three*nu2(i)*bb*b_i
652 cq1=c1*b
653 cq2=c1*db*2.0
654 c1=zep444*gama2(i)
655 dfn=an(i)*(cp1*dfnp*c1-fnp*cp2)
656 . + bn(i)*(cq1*dfnq*c1-fnq*cq2)
657 dfm=am(i)*(cp1*dfmp-fmp*cp2)
658 . + bm(i)*(cq1*dfmq-fmq*cq2)
659 dfnm=anm(i)*(cp1*dfnmp-fnmp*cp2)
660 . + bnm(i)*(cq1*dfnmq-fnmq*cq2)
661 df =(dfn+dfm+s*dfnm)*
662 . (e(i)*half-dr*h(i))/yld_i-2.*h(i)*yld_i
663
664
665C DEBUG
666C IF(I==21) THEN
667C WRITE(12,*) 'DFN,DFM,DFNM'
668C WRITE(12,*) DFN,DFM,DFNM
669C ERR=ABS(F/(DF*DPLA_I))
670C WRITE(12,*) 'N,ERR,DPLA_I,F,DF'
671C WRITE(12,'(I5,F8.5,3E16.6)') N,ERR,DPLA_I,F,DF
672C ENDIF
673C
674 dpla_j(i)=max(zero,dpla_i-f/df)
675 ENDDO
676 ENDDO
677C------------------------------------------
678C CONTRAINTES PLASTIQUEMENT ADMISSIBLES
679C------------------------------------------
680#include "vectorize.inc"
681 DO j=1,nindx
682 i=index(j)
683 pla(i) = pla(i) + dpla_j(i)
684 dpla_i=dpla_j(i)
685 yld_i =yld(i)+h(i)*dpla_i
686 dr =half*e(i)*dpla_i/yld_i
687 xp =dr*nu1(i)
688 xq =three*dr*nu2(i)
689 xpg =xp*zep444*gama2(i)
690 xqg =xq*zep444*gama2(i)
691 c1=one+qtier(i)
692 a=one+c1*xp+twop444*gama2(i)*xp*xp
693 b=one+c1*xq+twop444*gama2(i)*xq*xq
694 a_i=one/a
695 b_i=one/b
696 aa=a_i*a_i
697 bb=b_i*b_i
698 fnmp=one-onep222*gama2(i)*xp*xp
699 fnmq=one-onep222*gama2(i)*xq*xq
700 fnm=aa*fnmp*anm(i)+bb*fnmq*bnm(i)
701 IF (fnm<zero) THEN
702 s=-one
703 ELSE
704 s=one
705 ENDIF
706 pn=(one+qtier(i)*xp)*a_i
707 p_m=(one+xp)*a_i
708 pnm1=-sqr4_3*gama(i)*s*xp*a_i
709 pnm2=pnm1*one_over_12
710 qn=(one+qtier(i)*xq)*b_i
711 qm=(one+xq)*b_i
712 qnm1=-sqr4_3*xq*gama(i)*s*b_i
713 qnm2=qnm1*one_over_12
714 sn1=s1(i)*pn+sm1(i)*pnm1
715 sn2=s2(i)*qn+sm2(i)*qnm1
716 s3=signxy(i)*qn+momnxy(i)*qnm1
717 m1=sm1(i)*p_m+s1(i)*pnm2
718 m2=sm2(i)*qm+s2(i)*qnm2
719 momnxy(i)=signxy(i)*qnm2+momnxy(i)*qm
720 signxx(i)=sn1+sn2
721 signyy(i)=sn1-sn2
722 signxy(i)=s3
723 momnxx(i)=m1+m2
724 momnyy(i)=m1-m2
725 dezz = - nu3(i)*dr*sn1*2./e(i)
726 thk(i) = thk(i) + dezz*thk(i)*off(i)
727 ENDDO
728 ELSE
729C-------------------------
730C ITERATIVE
731C-------------------------
732C-------------------------
733C CRITERE DE VON MISES
734C-------------------------
735 DO i=1,nel
736C-------------------------------------------------------------------------
737C GAMA (L'INVERSE DE GAMA DANS LA FORMULE)
738C C GAMA (L'INVERSE DE GAMA DANS LA FORMULE) A ETE PRIS A 2/3 PAR DEFAUT
739C-------------------------------------------------------------------------
740 c1 = pla(i)*e(i)
741 ccc=exp(-twop6667*c1/yld(i))
742 gama(i) = two/(three-ccc)
743 gama2(i) = gama(i)*gama(i)
744 cm(i) = thirty6*gama2(i)
745 cnm(i) = threep4641*gama(i)
746 qtier(i) = three*gama2(i)
747 h(i) = max(zero,h(i))
748 s1(i) = (signxx(i)+signyy(i))*half
749 s2(i) = (signxx(i)-signyy(i))*half
750 s3 = signxy(i)
751 sm1(i) = (momnxx(i)+momnyy(i))*half
752 sm2(i) = (momnxx(i)-momnyy(i))*half
753 sm3 = momnxy(i)
754 an(i) = s1(i)*s1(i)
755 bn(i) = three*(s2(i)*s2(i)+s3*s3)
756 nvm(i) = an(i)+bn(i)
757 am(i) = sm1(i)*sm1(i)*cm(i)
758 bm(i) = three*(sm2(i)*sm2(i)+sm3*sm3)*cm(i)
759 mvm(i) = am(i)+bm(i)
760 anm(i) = s1(i)*sm1(i)*cnm(i)
761 bnm(i) = three*(s2(i)*sm2(i)+s3*sm3)*cnm(i)
762 nmvm(i) = anm(i)+bnm(i)
763 svm(i) = sqrt(nvm(i)+mvm(i)+abs(nmvm(i)))
764 dezz = -(depsxx(i)+depsyy(i))*nnu1
765 thk(i) = thk(i) +thk(i)* dezz*off(i)
766 ENDDO
767C-------------------------
768C GATHER PLASTIC FLOW
769C-------------------------
770 nindx=0
771 DO i=1,nel
772 IF(svm(i)>yld(i).AND.off(i)==one) THEN
773 nindx=nindx+1
774 index(nindx)=i
775 ENDIF
776 ENDDO
777 IF(nindx==0) RETURN
778C---------------------------
779C DEP EN CONTRAINTE PLANE
780C---------------------------
781 DO j=1,nindx
782 i=index(j)
783 nu=uparam(iadbufv+6)
784 nu1(i) = half*(one+nu)
785 nu2(i) = three_half*(one -nu)
786 nu3(i) = one -nnu1
787 num1(i) = one +qtier(i)
788 num2(i) = fivep5*gama2(i)
789 lfn(i)=num2(i)
790 qfn(i)=sixteenp5*gama2(i)*gama2(i)
791 qfnm(i)=-num2(i)
792 dpla_j(i)=(svm(i)-yld(i))/(g3(i)*qtier(i)+h(i))
793 etse(i)= h(i)/(h(i)+e(i))
794 ENDDO
795C-------------------------------
796C TIENT COMPTE DU COUPLAGE
797C-------------------------------
798 DO n=1,nmax
799#include "vectorize.inc"
800 DO j=1,nindx
801 i=index(j)
802 dpla_i=dpla_j(i)
803 yld_i =yld(i)+h(i)*dpla_i
804 dr =a1(i)*dpla_i/yld_i
805 xp =dr*nu1(i)
806 xq =dr*nu2(i)
807 da=num1(i)+num2(i)*xp
808 db=num1(i)+num2(i)*xq
809 a=one+(da+num1(i))*xp*half
810 b=one+(db+num1(i))*xq*half
811 a_i=one/a
812 b_i=one/b
813 aa=a_i*a_i
814 bb=b_i*b_i
815 dfnp=lfn(i)+qfn(i)*xp
816 dfnq=lfn(i)+qfn(i)*xq
817 dfmp=onep8333*(xp+one)
818 dfmq=onep8333*(xq+one)
819 dfnmp=qfnm(i)*xp
820 dfnmq=qfnm(i)*xq
821 xp = half*xp
822 xq = half*xq
823 fnp=one+(dfnp+lfn(i))*xp
824 fnq=one+(dfnp+lfn(i))*xq
825 fmp=one+(dfmp+onep8333)*xp
826 fmq=one+(dfmq+onep8333)*xq
827 fnmp=one+dfnmp*xp
828 fnmq=one+dfnmq*xq
829 fnm=aa*fnmp*anm(i)+bb*fnmq*bnm(i)
830 IF (fnm<zero) THEN
831 s=-one
832 ELSE
833 s=one
834 ENDIF
835 cp1 =(fnp*an(i)+s*fnmp*anm(i)+fmp*am(i))*aa
836 cq1 =(fnq*bn(i)+s*fnmq*bnm(i)+fmq*bm(i))*bb
837 cp2 =(dfnp*an(i)+s*dfnmp*anm(i)+dfmp*am(i))*aa
838 cq2 =(dfnq*bn(i)+s*dfnmq*bnm(i)+dfmq*bm(i))*bb
839 xpg =two*nu1(i)*da*a_i
840 xqg =two*nu2(i)*db*b_i
841 f =cp1 +cq1-yld_i*yld_i
842 df =(cp2*nu1(i)+cq2*nu2(i)-cp1*xpg-cq1*xqg)*
843 . (a1(i)-dr*h(i))/yld_i- two*h(i)*yld_i
844 dpla_j(i)=max(zero,dpla_i-f/df)
845 ENDDO
846 ENDDO
847C------------------------------------------
848C CONTRAINTES PLASTIQUEMENT ADMISSIBLES
849C------------------------------------------
850#include "vectorize.inc"
851 DO j=1,nindx
852 i=index(j)
853 pla(i) = pla(i) + dpla_j(i)
854 dpla_i=dpla_j(i)
855 yld_i =yld(i)+h(i)*dpla_i
856 dr =a1(i)*dpla_i/yld_i
857 xp =dr*nu1(i)
858 xq =dr*nu2(i)
859 xpg =xp*xp
860 xqg =xq*xq
861 a=one + num1(i)*xp+num2(i)*xpg
862 b=one + num1(i)*xq+num2(i)*xqg
863 a_i=one/a
864 b_i=one/b
865 aa=a_i*a_i
866 bb=b_i*b_i
867 fnmp=one + qfnm(i)*xpg
868 fnmq=one + qfnm(i)*xqg
869 fnm=aa*fnmp*anm(i)+bb*fnmq*bnm(i)
870 IF (fnm<zero) THEN
871 s=-onep732*gama(i)
872 ELSE
873 s=onep732*gama(i)
874 ENDIF
875 qn=one+qtier(i)*xq
876 qnm1=xq*s
877 qnm2=qnm1*one_over_12
878 sn1=(s1(i)*(one + qtier(i)*xp)-sm1(i)*s*xp)*a_i
879 sn2=(s2(i)*qn-sm2(i)*qnm1)*b_i
880 s3=(signxy(i)*qn-momnxy(i)*qnm1)*b_i
881 m1=(sm1(i)*(one + xp)-s1(i)*s*xp*one_over_12)*a_i
882 m2=(sm2(i)*(one + xq)-s2(i)*qnm2)*b_i
883 momnxy(i)=(momnxy(i)*(1.+xq)-signxy(i)*qnm2)*b_i
884 signxx(i)=sn1+sn2
885 signyy(i)=sn1-sn2
886 signxy(i)=s3
887 momnxx(i)=m1+m2
888 momnyy(i)=m1-m2
889 dezz = - nu3(i)*dr*sn1/e(i)
890 thk(i) = thk(i) + dezz*thk(i)*off(i)
891 ENDDO
892 ENDIF
893C
894 DO i=1,nel
895 IF(pla(i)>epsmax.AND.off(i)==one) THEN
896 off(i)=four_over_5
897 ENDIF
898 ENDDO
899C
900 RETURN
901 END
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
inter_rat
subroutine inter_rat(x0, x1, x2, x3, y0, y1, y2, y3, x, y, yp, i, n)
Definition sigeps60c.F:786
sigeps60g
subroutine sigeps60g(nel, nuvar, npf, tf, time, uparam, rho0, area, eint, thk0, epspxx, epspyy, epspxy, epspyz, epspzx, depsxx, depsyy, depsxy, depsyz, depszx, depbxx, depbyy, depbxy, epsxx, epsyy, epsxy, epsyz, epszx, sigoxx, sigoyy, sigoxy, sigoyz, sigozx, momoxx, momoyy, momoxy, signxx, signyy, signxy, signyz, signzx, momnxx, momnyy, momnxy, sigvxx, sigvyy, sigvxy, sigvyz, sigvzx, soundsp, viscmax, thk, pla, uvar, off, ngl, ipm, mat, etse, gs, yld, epsp, israte, ipla, shf)
Definition sigeps60g.F:49
vinter
subroutine vinter(tf, iad, ipos, ilen, nel, x, dydx, y)
Definition vinter.F:72