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