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