OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps36g.F File Reference
#include "implicit_f.inc"
#include "mvsiz_p.inc"
#include "param_c.inc"
#include "scr03_c.inc"
#include "scr05_c.inc"
#include "vectorize.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine sigeps36g (nel, 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, off, ngl, ipm, mat, etse, gs, yld, epsp, israte, ipla, shf, nvartmp, vartmp)

Function/Subroutine Documentation

◆ sigeps36g()

subroutine sigeps36g ( integer nel,
integer, dimension(*) 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,
off,
integer, dimension(nel) ngl,
integer, dimension(npropmi,*) ipm,
integer, dimension(nel) mat,
etse,
gs,
yld,
epsp,
integer israte,
integer ipla,
shf,
integer nvartmp,
integer, dimension(nel,nvartmp), intent(inout) vartmp )

Definition at line 31 of file sigeps36g.F.

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