Télécharger behav3.eso

Retour à la liste

Numérotation des lignes :

  1. C BEHAV3 SOURCE CB215821 16/04/21 21:15:19 8920
  2. SUBROUTINE BEHAV3(SIGR,DSTRN,DEPS1,DEPS2,IPLA,
  3. & SIG1,SIG2,IFIS,SIGF,DSIGT,NSTRS,IFOUB,DEP,SIGRV,SIGP,
  4. & BETJEF,VISCO,NECH0,NECH1)
  5. C
  6. C ==================================================================
  7. C
  8. C MODELE DE PLASTICITE EN TRAITEMENT PRE ET POST PIC
  9. C Deux criteres de traction compression:
  10. C Rankine et Von Mises
  11. C
  12. C ==================================================================
  13. C CE SOUS-PROGRAMME EST APPELE DANS "BONE".
  14. C
  15. IMPLICIT INTEGER(I-N)
  16. IMPLICIT REAL*8(A-H,O-Z)
  17. DIMENSION DFSIG1(4),DFSIG2(4),SIGF(4),DSIGT(4),VEC1(4),VEC2(4)
  18. DIMENSION V1(4),AC1(4),AC2(4),SIGE(4),SIGP(4)
  19. DIMENSION DLAM1(2),DLAM0(2),SIGR(4),IS1(2)
  20. DIMENSION FCRI0(2),FCRI1(2),DEPSI(4),DSTRN(4),SIGRV(4)
  21. DIMENSION D2FSI1(4,4),D2FSI2(4,4),DEP(4,4),P1(4,4),P2(4,4)
  22. DIMENSION D(4,4),A(4,4),AI(4,4),DJAC0(4,4),DJI(4,4),AH(4,4)
  23. C
  24. *
  25. SEGMENT BETJEF
  26. REAL*8 AA,BETA,COLI,PALF,YOUN,XNU,GFC,GFT,CAR,ETA,TDEF,
  27. & TCON,DPSTF1,DPSTF2,TETA,PDT,TP0
  28. INTEGER ICT,ICC,IMOD,IVIS,ITR,
  29. & ISIM,IBB,IGAU,IZON
  30. ENDSEGMENT
  31. SEGMENT VISCO
  32. REAL*8 DPSTV1,DPSTV2,SIGV1,SIGV2,ENDV
  33. ENDSEGMENT
  34. SEGMENT NECH0
  35. REAL*8 DT,DC,ALFG,S0,ENDO
  36. ENDSEGMENT
  37. SEGMENT NECH1
  38. REAL*8 ENDL
  39. ENDSEGMENT
  40. *
  41. *
  42. * COMMON /DBETJEF/AA,BETA,COLI,PALF,YOUN,XNU,GFC,GFT,CAR,ETA,TDEF,
  43. * & TCON,DPSTF1,DPSTF2,TETA,PDT,ICT,ICC,IMOD,IVIS,ITR,
  44. * & ISIM,IBB,IGAU,IZON
  45. C
  46. MNSTRS=NSTRS
  47. C
  48. IAPEX=0
  49. Ft=PALF*COLI
  50. Rb=COLI
  51. PRB=1.D-10
  52. PRB2=1.D-6
  53. ITER=1
  54. IEC1=0
  55. IEC2=0
  56. IBROY=0
  57. ITANG=0
  58. CRIMAX=0.D0
  59. SEQ = 0.D0
  60. SEQ1 = 0.D0
  61. SEQ2 = 0.D0
  62. CALL ZERO(SIGE,4,1)
  63. CALL ZERO(V1,4,1)
  64. CALL ZERO(P1,4,4)
  65. CALL ZERO(P2,4,4)
  66. CALL ZERO(D,4,4)
  67. C
  68. DO 10 I=1,NSTRS
  69. SIGE(I)=SIGF(I)
  70. 10 CONTINUE
  71. C
  72. IF (IMOD.EQ.1.OR.IMOD.EQ.3) THEN
  73. AD=YOUN/(1.D0-XNU*XNU)
  74. D(1,1)=AD
  75. D(2,2)=D(1,1)
  76. D(3,3)=AD*(1.D0-XNU)/2.D0
  77. D(1,2)=AD*XNU
  78. D(2,1)=D(1,2)
  79. ENDIF
  80. C
  81. IF (IMOD.EQ.2.OR.IMOD.EQ.4) THEN
  82. ADD=YOUN/((1.D0+XNU)*(1.D0-2.D0*XNU))
  83. D(1,1)=ADD*(1.D0-XNU)
  84. D(2,2)=D(1,1)
  85. D(3,3)=D(1,1)
  86. D(1,2)=ADD*XNU
  87. D(2,1)=D(1,2)
  88. D(1,3)=D(1,2)
  89. D(2,3)=D(1,2)
  90. D(3,1)=D(1,2)
  91. D(3,2)=D(1,2)
  92. D(4,4)=0.5*ADD*(1.D0-2.D0*XNU)
  93. ENDIF
  94. C
  95. C ************ Le point est fissure pour 1ere fois ***********
  96. C
  97. IF (IFIS.EQ.0) THEN
  98. CALL PRINC(SIGF,V1,NSTRS)
  99. DEPS1=0.D0
  100. SIG1=Ft
  101. IF (V1(1).GT.Ft) THEN
  102. IFIS=1
  103. ENDIF
  104. ENDIF
  105. IF (IPLA.EQ.0) THEN
  106. IF(IMOD.EQ.1.OR.IMOD.EQ.2) THEN
  107. CALL CRIVON(SIGF,SEQ,NSTRS,BETJEF)
  108. ENDIF
  109. IF(IMOD.EQ.3.OR.IMOD.EQ.4) THEN
  110. CALL CRIDRU(SIGF,SEQ,NSTRS,BETJEF)
  111. ENDIF
  112. DEPS2=0.D0
  113. SIG2=Rb*AA
  114. IF (SEQ.GT.SIG2) THEN
  115. IPLA=1
  116. ENDIF
  117. ENDIF
  118. C
  119. 12 CONTINUE
  120. C
  121. CALL PRINC(SIGE,V1,NSTRS)
  122. TETA=V1(4)
  123. IF(IMOD.EQ.1.OR.IMOD.EQ.2) THEN
  124. CALL CRIVON(SIGE,SEQ,NSTRS,BETJEF)
  125. ENDIF
  126. IF(IMOD.EQ.3.OR.IMOD.EQ.4) THEN
  127. CALL CRIDRU(SIGE,SEQ,NSTRS,BETJEF)
  128. ENDIF
  129. FCRI0(1)=V1(1)-SIG1
  130. FCRI0(2)=SEQ-SIG2
  131. IF (ABS(FCRI0(1)).GT.ABS(FCRI0(2))) THEN
  132. CRIMAX=ABS(100.D0*FCRI0(1))
  133. ELSE
  134. CRIMAX=ABS(100.D0*FCRI0(2))
  135. ENDIF
  136. IF (FCRI0(1).LT.0.D0.AND.FCRI0(2).LT.0.D0) THEN
  137. IEC1=1
  138. IEC2=1
  139. ENDIF
  140. IF (FCRI0(1).LT.0.D0.AND.FCRI0(2).GE.0.D0) THEN
  141. IEC1=1
  142. IEC2=0
  143. ENDIF
  144. IF (FCRI0(1).GE.0.D0.AND.FCRI0(2).LT.0.D0) THEN
  145. IEC1=0
  146. IEC2=1
  147. ENDIF
  148. C
  149. 15 CONTINUE
  150. C
  151. DO 16 I=1,NSTRS
  152. SIGF(I)=SIGE(I)
  153. 16 CONTINUE
  154. C
  155. IF (IEC1.EQ.1.AND.IEC2.EQ.1) THEN
  156. DO 17 I=1,NSTRS
  157. DO 17 J=1,NSTRS
  158. DEP(I,J)=D(I,J)
  159. 17 CONTINUE
  160. GOTO 100
  161. ENDIF
  162. C
  163. IF (IEC1.EQ.1.AND.IEC2.EQ.0) THEN
  164. CALL BEHAV2(SIGR,DSTRN,DEPS1,DEPS2,IPLA
  165. & ,SIG1,SIG2,IFIS,SIGF,DSIGT,NSTRS,IFOUB,DEP,SIGRV,SIGP,
  166. & BETJEF,VISCO,NECH0,NECH1)
  167. GOTO 100
  168. ENDIF
  169. C
  170. IF (IEC1.EQ.0.AND.IEC2.EQ.1) THEN
  171. CALL BEHAV1(SIGR,DSTRN,DEPS1,DEPS2,IPLA
  172. & ,SIG1,SIG2,IFIS,SIGF,DSIGT,NSTRS,IFOUB,DEP,SIGRV,SIGP,
  173. & BETJEF,VISCO,NECH0,NECH1)
  174. GOTO 100
  175. ENDIF
  176. C
  177. DK1=DEPS1
  178. DK2=DEPS2
  179. DLAM0(1)=0.D0
  180. DLAM0(2)=0.D0
  181. C
  182. 18 CONTINUE
  183. C
  184. C ************ Preparation des criteres ********************
  185. C
  186. C ------------------ Von Mises -----------------------------
  187. IF (IMOD.EQ.1.OR.IMOD.EQ.2) THEN
  188. A1=0.D0
  189. A2=1.D0
  190. ENDIF
  191. C ----------------- Drucker Prager -------------------------
  192. IF (IMOD.EQ.3.OR.IMOD.EQ.4) THEN
  193. A1=(BETA-1.D0)/(2*BETA-1.D0)
  194. A2=BETA/(2*BETA-1.D0)
  195. ENDIF
  196. C ----------------------------------------------------------
  197. C
  198. C ---------------- direction de traction ------------------
  199. C
  200. C--------------- Determination du parametre d'ecrouissage ----------
  201. C
  202. IF(IVIS.LE.2) THEN
  203. CALL UNICOU(DK1,PAECT,1,SEQ1,BETJEF)
  204. ELSE
  205. CALL UNICO1(DK1,PAECT,1,SEQ1,BETJEF,NECH0,NECH1)
  206. ENDIF
  207. C
  208. C------------------------------------------------------------
  209. PI=4.D0*ATAN(1.D0)
  210. PHIC=V1(4)*(PI/180.D0)
  211. COSA=COS(PHIC)
  212. SINA=SIN(PHIC)
  213. C
  214. IF (IMOD.EQ.1.OR.IMOD.EQ.3) THEN
  215. DFSIG1(1)=COSA*COSA
  216. DFSIG1(2)=SINA*SINA
  217. DFSIG1(3)=2.D0*SINA*COSA
  218. DO 19 I=1,NSTRS
  219. DO 19 J=1,NSTRS
  220. P1(I,J)=0.D0
  221. 19 CONTINUE
  222. P1(1,1)=1.D0/2.D0
  223. P1(1,2)=-1.D0/2.D0
  224. P1(2,1)=-1.D0/2.D0
  225. P1(2,2)=1.D0/2.D0
  226. P1(3,3)=2.D0
  227. ENDIF
  228. C
  229. IF (IMOD.EQ.2.OR.IMOD.EQ.4) THEN
  230. DFSIG1(1)=COSA*COSA
  231. DFSIG1(2)=SINA*SINA
  232. DFSIG1(3)=0.D0
  233. DFSIG1(4)=2.D0*SINA*COSA
  234. DO 21 I=1,NSTRS
  235. DO 21 J=1,NSTRS
  236. P1(I,J)=0.D0
  237. 21 CONTINUE
  238. P1(1,1)=1.D0/2.D0
  239. P1(1,2)=-1.D0/2.D0
  240. P1(2,1)=-1.D0/2.D0
  241. P1(2,2)=1.D0/2.D0
  242. P1(4,4)=2.D0
  243. ENDIF
  244. C
  245. DO 20 I=1,NSTRS
  246. AC1(I)=0.D0
  247. DO 20 J=1,NSTRS
  248. AC1(I)=AC1(I)+D(I,J)*DFSIG1(J)
  249. 20 CONTINUE
  250. C
  251. C ---------------- direction de compression ----------------
  252. C
  253. IF(IVIS.LE.2) THEN
  254. CALL UNICOU(DK2,PAEC,2,SEQ2,BETJEF)
  255. ELSE
  256. CALL UNICO1(DK2,PAEC,2,SEQ2,BETJEF,NECH0,NECH1)
  257. ENDIF
  258. C
  259. IF (IMOD.EQ.1.OR.IMOD.EQ.3) THEN
  260. SX=SIGF(1)
  261. SY=SIGF(2)
  262. SXY=SIGF(3)
  263. R1=(SX-SY)*(SX-SY)+SX*SX+SY*SY+6*SXY*SXY
  264. R=SQRT(0.5*R1)
  265. IF (ABS(R).LE.1.D-10) THEN
  266. WRITE(*,*)'INDETERMINATION de DFSIG ds BEHAV3'
  267. STOP
  268. ENDIF
  269. DFSIG2(1)=(2.D0*SX-SY)/2.D0/R/A2+A1/A2
  270. DFSIG2(2)=(2.D0*SY-SX)/2.D0/R/A2+A1/A2
  271. DFSIG2(3)=3.D0*SXY/R/A2
  272. P2(1,1)=2.D0/A2/A2
  273. P2(1,2)=-1.D0/A2/A2
  274. P2(1,3)=0.D0/A2/A2
  275. P2(2,1)=-1.D0/A2/A2
  276. P2(2,2)=2.D0/A2/A2
  277. P2(2,3)=0.D0/A2/A2
  278. P2(3,1)=0.D0/A2/A2
  279. P2(3,2)=0.D0/A2/A2
  280. P2(3,3)=6.D0/A2/A2
  281. ENDIF
  282. C
  283. IF (IMOD.EQ.2.OR.IMOD.EQ.4) THEN
  284. SX=SIGF(1)
  285. SY=SIGF(2)
  286. SZ=SIGF(3)
  287. SXY=SIGF(4)
  288. TX=(SX-SY)*(SX-SY)
  289. TY=(SX-SZ)*(SX-SZ)
  290. TZ=(SY-SZ)*(SY-SZ)
  291. R1=TX+TY+TZ+6*SXY*SXY
  292. R=SQRT(0.5*R1)
  293. IF (ABS(R).LE.1.D-10) THEN
  294. WRITE(*,*)'INDETERMINATION de DFSIG ds BEHAV3'
  295. STOP
  296. ENDIF
  297. DFSIG2(1)=(2.D0*SX-SY-SZ)/2.D0/R/A2+A1/A2
  298. DFSIG2(2)=(2.D0*SY-SX-SZ)/2.D0/R/A2+A1/A2
  299. DFSIG2(3)=(2.D0*SZ-SX-SY)/2.D0/R/A2+A1/A2
  300. DFSIG2(4)=3.D0*SXY/R/A2
  301. P2(1,1)=2.D0/A2/A2
  302. P2(1,2)=-1.D0/A2/A2
  303. P2(1,3)=-1.D0/A2/A2
  304. P2(1,4)=0.D0/A2/A2
  305. P2(2,1)=-1.D0/A2/A2
  306. P2(2,2)=2.D0/A2/A2
  307. P2(2,3)=-1.D0/A2/A2
  308. P2(2,4)=0.D0/A2/A2
  309. P2(3,1)=-1.D0/A2/A2
  310. P2(3,2)=-1.D0/A2/A2
  311. P2(3,3)=2.D0/A2/A2
  312. P2(3,4)=0.D0/A2/A2
  313. P2(4,1)=0.D0/A2/A2
  314. P2(4,2)=0.D0/A2/A2
  315. P2(4,3)=0.D0/A2/A2
  316. P2(4,4)=6.D0/A2/A2
  317. ENDIF
  318. C
  319. DO 25 I=1,NSTRS
  320. AC2(I)=0.D0
  321. DO 25 J=1,NSTRS
  322. AC2(I)=AC2(I)+D(I,J)*DFSIG2(J)
  323. 25 CONTINUE
  324. C
  325. C ************ Preparation du jacobien ********************
  326. C
  327. F1DF1=0.D0
  328. DO 30 J=1,NSTRS
  329. F1DF1=F1DF1+DFSIG1(J)*AC1(J)
  330. 30 CONTINUE
  331. F1DF2=0.D0
  332. DO 31 J=1,NSTRS
  333. F1DF2=F1DF2+DFSIG1(J)*AC2(J)
  334. 31 CONTINUE
  335. F2DF2=0.D0
  336. DO 32 J=1,NSTRS
  337. F2DF2=F2DF2+DFSIG2(J)*AC2(J)
  338. 32 CONTINUE
  339. F2DF1=0.D0
  340. DO 33 J=1,NSTRS
  341. F2DF1=F2DF1+DFSIG2(J)*AC1(J)
  342. 33 CONTINUE
  343. C
  344. DJAC0(1,1)=-(F1DF1+PAECT)
  345. DJAC0(1,2)=-(F1DF2)
  346. DJAC0(2,2)=-(F2DF2+PAEC)
  347. DJAC0(2,1)=-(F2DF1)
  348. C
  349. DO 43 I=1,2
  350. DO 43 J=1,2
  351. DJI(I,J)=DJAC0(I,J)
  352. 43 CONTINUE
  353. C
  354. CALL INVMA2(DJI,2,ISING)
  355. IF (ISING.EQ.1) THEN
  356. WRITE(*,*)'MATRICE DJI singuliere ds BEHAV3'
  357. ENDIF
  358. C
  359. C ************ Debut des iterations internes *******************
  360. C
  361. 40 CONTINUE
  362. C
  363. C *************** Determination de DK et DLAM ******************
  364. C
  365. C
  366. C
  367. DLAM1(1)=DLAM0(1)-DJI(1,1)*FCRI0(1)-DJI(1,2)*FCRI0(2)
  368. DLAM1(2)=DLAM0(2)-DJI(2,1)*FCRI0(1)-DJI(2,2)*FCRI0(2)
  369. C
  370. IF (DLAM1(1).LE.0.D0) IEC1=1
  371. IF (DLAM1(2).LE.0.D0) IEC2=1
  372. IF (IEC1.EQ.1.OR.IEC2.EQ.1) THEN
  373. C WRITE(*,*)'Dans BEHAV3, DLAMDA1 est negatif:',DLAM1(1)
  374. C WRITE(*,*)'Dans BEHAV3, DLAMDA2 est negatif:',DLAM1(2)
  375. C WRITE(*,*)'A l iteration :',ITER
  376. GOTO 15
  377. ENDIF
  378. C
  379. DK1=DEPS1+DLAM1(1)
  380. DK2=DEPS2+DLAM1(2)
  381. C
  382. IF(IVIS.LE.2) THEN
  383. CALL UNICOU(DK1,PAECT,1,SEQ1,BETJEF)
  384. CALL UNICOU(DK2,PAEC,2,SEQ2,BETJEF)
  385. ELSE
  386. CALL UNICO1(DK1,PAECT,1,SEQ1,BETJEF,NECH0,NECH1)
  387. CALL UNICO1(DK2,PAEC,2,SEQ2,BETJEF,NECH0,NECH1)
  388. ENDIF
  389. C
  390. C ************** Determination de DPHI *********************
  391. C
  392. IF (IMOD.EQ.1) THEN
  393. AD1=YOUN/(1.D0-XNU)
  394. coef=2.D0*SEQ2/(2.D0*SEQ2+DLAM1(2)*AD1)
  395. DPHI1=SEQ1-0.5*coef*((SIGE(1)+SIGE(2))-AD1*DLAM1(1))
  396. TO=SIGF(3)
  397. DPHI2=SEQ2
  398. ENDIF
  399. IF (IMOD.EQ.2) THEN
  400. DG=YOUN/2.D0/(1.D0+XNU)
  401. AD1=YOUN/(1.D0+XNU)/(1.D0-2.D0*XNU)
  402. coef1=(SEQ2+2.D0*DLAM1(2)*DG)/(SEQ2+3.D0*DLAM1(2)*DG)
  403. coef2=(2.D0*DLAM1(2)*DG)/(SEQ2+3.D0*DLAM1(2)*DG)
  404. CAS1=(SIGE(1)+SIGE(2)-0.5*DLAM1(1)*2.D0*AD1)
  405. CAS2=SIGE(3)-DLAM1(1)*XNU*AD1
  406. DPHI1=SEQ1-0.5*(coef1*CAS1+coef2*CAS2)
  407. TO=SIGF(4)
  408. DPHI2=SEQ2
  409. ENDIF
  410. IF (IMOD.EQ.3) THEN
  411. DPHI2=SEQ2-A1/A2*(SIGF(1)+SIGF(2))
  412. AD1=YOUN/(1.D0-XNU)
  413. coef=2.D0*DPHI2/(2.D0*DPHI2+DLAM1(2)*AD1)
  414. SOM=A1/A2*DLAM1(2)*2.D0*AD1
  415. DPHI1=SEQ1-0.5*coef*((SIGE(1)+SIGE(2))-AD1*DLAM1(1)-SOM)
  416. TO=SIGF(3)
  417. ENDIF
  418. IF (IMOD.EQ.4) THEN
  419. SOM=SIGE(1)+SIGE(2)+SIGE(3)
  420. DPHI2=A1*SOM/A2
  421. DPHI2=DPHI2-A1*DLAM1(1)*ADD*(1.D0+XNU)/A2
  422. DPHI2=DPHI2-A1*A1*DLAM1(2)*3.D0*ADD*(1.D0+XNU)/A2/A2
  423. DPHI2=SEQ2-DPHI2
  424. DG=YOUN/2.D0/(1.D0+XNU)/A2/A2
  425. AD1=YOUN/(1.D0+XNU)/(1.D0-2.D0*XNU)
  426. AD2=YOUN/(1.D0-2.D0*XNU)
  427. coef1=(DPHI2+2.D0*DLAM1(2)*DG)/(DPHI2+3.D0*DLAM1(2)*DG)
  428. coef2=(2.D0*DLAM1(2)*DG)/(DPHI2+3.D0*DLAM1(2)*DG)
  429. CAS1=(SIGE(1)+SIGE(2)-DLAM1(1)*AD1-2.D0*A1/A2*DLAM1(2)*AD2)
  430. CAS2=SIGE(3)-DLAM1(1)*XNU*AD1-A1/A2*DLAM1(2)*AD2
  431. DPHI1=SEQ1-0.5*(coef1*CAS1+coef2*CAS2)
  432. TO=SIGF(4)
  433. ENDIF
  434. C IF (DPHI1.LT.0.D0) THEN
  435. C WRITE(*,*)'ATTENTION DPHI1 NEGATIF'
  436. C WRITE(*,*)'DPHI1=',DPHI1
  437. C ENDIF
  438. C IF (DPHI2.LT.0.D0) THEN
  439. C WRITE(*,*)'ATTENTION DPHI2 ds BEHAV3 NEGATIF'
  440. C WRITE(*,*)'DPHI2=',DPHI2
  441. C ENDIF
  442. C
  443. C--------------- Cas de l'apex -----------------------------
  444. C
  445. IF (ABS(DPHI1).LE.10E-10.AND.
  446. *ABS(DPHI2).GT.10E-10) THEN
  447. IAPEX=1
  448. C WRITE(*,*)'IAPEX ds BEHAV3 =',IAPEX
  449. C WRITE(*,*)'Dans l element',IBB
  450. C WRITE(*,*)'et au point d intégration',IGAU
  451. DO 50 I=1,NSTRS
  452. DO 50 J=1,NSTRS
  453. AI(I,J)=0.D0
  454. 50 CONTINUE
  455. AI(1,1)=0.5
  456. AI(1,2)=0.5
  457. AI(2,1)=AI(1,2)
  458. AI(2,2)=AI(1,1)
  459. IF (IMOD.EQ.1.OR.IMOD.EQ.3) AI(3,3)=0.D0
  460. IF (IMOD.EQ.2.OR.IMOD.EQ.4) AI(3,3)=1.D0
  461. GOTO 75
  462. ENDIF
  463. C
  464. C---------- Cas du critere reduit a un point ----------------
  465. C
  466. IF (ABS(DPHI2).LE.10E-10) THEN
  467. IAPEX=2
  468. C WRITE(*,*)'IAPEX ds BEHAV3 =',IAPEX
  469. C WRITE(*,*)'Dans l element',IBB
  470. C WRITE(*,*)'et au point d intégration',IGAU
  471. DO 52 I=1,NSTRS
  472. DO 52 J=1,NSTRS
  473. AI(I,J)=0.D0
  474. 52 CONTINUE
  475. IF (IMOD.EQ.2.OR.IMOD.EQ.4) THEN
  476. DO 53 I=1,3
  477. DO 53 J=1,3
  478. AI(I,J)=1.D0/3.D0
  479. 53 CONTINUE
  480. ENDIF
  481. GOTO 75
  482. ENDIF
  483. C
  484. C ************** Mise a jour des contraintes ***************
  485. C
  486. C ---------------- calcul de la matrice A ------------------
  487. C
  488. DO 60 I=1,NSTRS
  489. DO 60 J=1,NSTRS
  490. A(I,J)=0.D0
  491. 60 CONTINUE
  492. C
  493. DG=YOUN/2.D0/(1.D0+XNU)
  494. C
  495. IF (IMOD.EQ.1.OR.IMOD.EQ.3) THEN
  496. DW1=AD*(2.D0-XNU)
  497. DW2=AD*(2.D0*XNU-1.D0)
  498. A(1,1)=1.D0+(DLAM1(1)*DG)/2.D0/DPHI1
  499. A(1,1)=A(1,1)+(DLAM1(2)*DW1)/2.D0/DPHI2/A2/A2
  500. A(2,2)=A(1,1)
  501. A(1,2)=-(DLAM1(1)*DG)/2.D0/DPHI1
  502. A(1,2)=A(1,2)+(DLAM1(2)*DW2)/2.D0/DPHI2/A2/A2
  503. A(2,1)=A(1,2)
  504. A(3,3)=(DLAM1(1)*DG)/DPHI1
  505. A(3,3)=A(3,3)+1.D0+(DLAM1(2)*3.D0*DG)/DPHI2/A2/A2
  506. ENDIF
  507. C
  508. IF (IMOD.EQ.2.OR.IMOD.EQ.4) THEN
  509. A(1,1)=1.D0+(DLAM1(1)*DG)/2.D0/DPHI1
  510. A(1,1)=A(1,1)+(DLAM1(2)*2.D0*DG)/DPHI2/A2/A2
  511. A(2,2)=A(1,1)
  512. A(3,3)=1.D0+(DLAM1(2)*2.D0*DG)/DPHI2/A2/A2
  513. A(1,2)=-(DLAM1(1)*DG)/2.D0/DPHI1
  514. A(1,2)=A(1,2)-(DLAM1(2)*DG)/DPHI2/A2/A2
  515. A(2,1)=A(1,2)
  516. A(1,3)=-(DLAM1(2)*DG)/DPHI2/A2/A2
  517. A(3,1)=A(1,3)
  518. A(3,2)=A(1,3)
  519. A(2,3)=A(1,3)
  520. A(4,4)=1.D0+(DLAM1(1)*DG)/DPHI1
  521. A(4,4)=A(4,4)+(DLAM1(2)*DG*3.D0)/DPHI2/A2/A2
  522. ENDIF
  523. C
  524. C -------------- invertion de la matrice A -----------------
  525. C
  526. DO 70 I=1,NSTRS
  527. DO 70 J=1,NSTRS
  528. AI(I,J)=A(I,J)
  529. 70 CONTINUE
  530. C
  531. CALL INVMA2(AI,NSTRS,ISING)
  532. IF (ISING.EQ.1) THEN
  533. WRITE(*,*)'MATRICE AI singuliere ds BEHAV3'
  534. ENDIF
  535. C
  536. C -------------- mise a jour des contraintes ------------
  537. C
  538. 75 CONTINUE
  539. VEC1(1)=AD1
  540. VEC1(2)=AD1
  541. VEC1(3)=0.D0
  542. VEC1(4)=0.D0
  543. IF(IMOD.EQ.3) THEN
  544. VEC2(1)=AD1*(1.D0+XNU)
  545. VEC2(2)=AD1*(1.D0+XNU)
  546. VEC2(3)=0.D0
  547. VEC2(4)=0.D0
  548. ENDIF
  549. IF(IMOD.EQ.4) THEN
  550. VEC2(1)=AD1*(1.D0+XNU)
  551. VEC2(2)=AD1*(1.D0+XNU)
  552. VEC2(3)=AD1*(1.D0+XNU)
  553. VEC2(4)=0.D0
  554. ENDIF
  555. IF(IMOD.EQ.1.OR.IMOD.EQ.2) THEN
  556. VEC2(1)=0.D0
  557. VEC2(2)=0.D0
  558. VEC2(3)=0.D0
  559. VEC2(4)=0.D0
  560. ENDIF
  561. C
  562. DO 80 I=1,NSTRS
  563. DEPSI(I)=SIGE(I)-0.5*DLAM1(1)*VEC1(I)
  564. DEPSI(I)=DEPSI(I)-A1*DLAM1(2)*VEC2(I)/A2
  565. 80 CONTINUE
  566. C
  567. DO 90 I=1,NSTRS
  568. SIGF(I)=0.0D+00
  569. DO 90 J=1,NSTRS
  570. SIGF(I)=SIGF(I)+AI(I,J)*DEPSI(J)
  571. 90 CONTINUE
  572. C
  573. C ******** Verification des criteres ****************
  574. C
  575. CALL PRINC(SIGF,V1,NSTRS)
  576. IF(IMOD.EQ.1.OR.IMOD.EQ.2) THEN
  577. CALL CRIVON(SIGF,SEQ,NSTRS,BETJEF)
  578. ENDIF
  579. IF(IMOD.EQ.3.OR.IMOD.EQ.4) THEN
  580. CALL CRIDRU(SIGF,SEQ,NSTRS,BETJEF)
  581. ENDIF
  582. FCRI1(1)=V1(1)-SEQ1
  583. IF (IAPEX.EQ.2) FCRI1(1)=0.D0
  584. FCRI1(2)=SEQ-SEQ2
  585. C
  586. IF (IBROY.EQ.0.AND.(ABS(FCRI1(1)).GE.CRIMAX.OR
  587. * .ABS(FCRI1(2)).GE.CRIMAX)) THEN
  588. C WRITE(*,*)'****************************************'
  589. C WRITE(*,*)'LE RESIDU DIVERGE AVEC BROYDEN'
  590. C WRITE(*,*)'on passe donc a la secante'
  591. C WRITE(*,*)'Dans l element',IBB
  592. C WRITE(*,*)'et au point d intégration',IGAU
  593. C WRITE(*,*)'CRIMAX=',CRIMAX
  594. C WRITE(*,*)'****************************************'
  595. ITER=ITR
  596. ENDIF
  597. C
  598. C ******* Compteur sur la methode de resolution ****
  599. C
  600. IF (IBROY.EQ.0.AND.ITER.EQ.ITR) THEN
  601. IBROY=1
  602. ITANG=1
  603. ITER=1
  604. IAPEX=0
  605. IEC1=0
  606. IEC2=0
  607. GOTO 12
  608. ENDIF
  609. C
  610. C ******* non convergence **************************
  611. C
  612. IF ((ABS(FCRI1(1)).GT.PRB.OR.ABS(FCRI1(2)).GT.PRB)
  613. *.AND.ITER.LT.ITR) THEN
  614. IF (IBROY.EQ.0) THEN
  615. CALL BROYDI(DJI,DLAM0,DLAM1,FCRI0,FCRI1)
  616. DLAM0(1)=DLAM1(1)
  617. DLAM0(2)=DLAM1(2)
  618. FCRI0(1)=FCRI1(1)
  619. FCRI0(2)=FCRI1(2)
  620. IF (ITER.GE.(ITR-1)) THEN
  621. C WRITE(*,*)'***********************'
  622. C WRITE(*,*)'BROYDEN n a pas aboutit'
  623. C WRITE(*,*)'ITER=',ITER
  624. C WRITE(*,*)'FCRIT=',FCRI0(1)
  625. C WRITE(*,*)'FCRIC=',FCRI0(2)
  626. C WRITE(*,*)'Dans l element',IBB
  627. C WRITE(*,*)'et au point d intégration',IGAU
  628. C WRITE(*,*)'***********************'
  629. ENDIF
  630. ITER=ITER+1
  631. GOTO 40
  632. ENDIF
  633. IF (IBROY.EQ.1.AND.ITANG.EQ.1) THEN
  634. DLAM0(1)=DLAM1(1)
  635. DLAM0(2)=DLAM1(2)
  636. FCRI0(1)=FCRI1(1)
  637. FCRI0(2)=FCRI1(2)
  638. CALL PRINC(SIGF,V1,NSTRS)
  639. IF (ITER.GE.(ITR-5)) THEN
  640. C WRITE(*,*)'ITER=',ITER
  641. C WRITE(*,*)'FCRIT=',FCRI0(1)
  642. C WRITE(*,*)'FCRIC=',FCRI0(2)
  643. ENDIF
  644. ITER=ITER+1
  645. GOTO 18
  646. ENDIF
  647. ENDIF
  648. IF (ITER.GE.ITR.AND.(ABS(FCRI1(1)).GT.PRB2.OR.
  649. * ABS(FCRI1(2)).GT.PRB2)) THEN
  650. WRITE(*,*)'NON CONVERGENCE INTERNE dans BEHAV3'
  651. WRITE(*,*)'Dans l element',IBB
  652. WRITE(*,*)'et au point d intégration',IGAU
  653. WRITE(*,*)'FCRIT=',FCRI0(1)
  654. WRITE(*,*)'FCRIC=',FCRI0(2)
  655. WRITE(*,*)'IPLA=',IPLA
  656. WRITE(*,*)'IFIS=',IFIS
  657. C STOP
  658. ENDIF
  659. C
  660. C **************** Fin des iterations internes *******************
  661. C
  662. C
  663. C ************* Calcul de la DEP ****************
  664. C
  665. IF (IAPEX.EQ.0) THEN
  666. CALL DERI2(D2FSI1,DPHI1,P1,SIGF,NSTRS,BETJEF)
  667. CALL DERI2(D2FSI2,DPHI2,P2,SIGF,NSTRS,BETJEF)
  668. CALL LADEP1(SIGF,DEP,PAECT,PAEC,NSTRS,IFOU,
  669. & DFSIG1,DFSIG2,D2FSI1,D2FSI2,DLAM1,BETJEF)
  670. ELSE IF (IAPEX.EQ.1) THEN
  671. DO 95 I=1,NSTRS
  672. DO 95 J=1,NSTRS
  673. AH(I,J)=0.D0
  674. DO 95 K=1,NSTRS
  675. AH(I,J)=AH(I,J)+AI(I,K)*D(K,J)
  676. 95 CONTINUE
  677. DFSIG1(1)=0.5*SQRT(2.D0)
  678. DFSIG1(2)=0.5*SQRT(2.D0)
  679. DFSIG1(3)=0.D0
  680. IF (IMOD.EQ.2) DFSIG1(4)=0.D0
  681. CALL CREDEP(AH,DFSIG1,PAEC,NSTRS,DEP,BETJEF)
  682. ELSE IF (IAPEX.EQ.2) THEN
  683. DO 96 I=1,NSTRS
  684. DO 96 J=1,NSTRS
  685. DEP(I,J)=0.D0
  686. DO 96 K=1,NSTRS
  687. DEP(I,J)=DEP(I,J)+AI(I,K)*D(K,J)
  688. 96 CONTINUE
  689. ENDIF
  690. IAPEX=0
  691. C DO 97 I=1,NSTRS
  692. C DO 97 J=1,NSTRS
  693. C DEP(I,J)=D(I,J)
  694. C 97 CONTINUE
  695. C
  696. C ***********************************************
  697. C
  698. DEPS1=DK1
  699. DEPS2=DK2
  700. SIG1=SEQ1
  701. SIG2=SEQ2
  702. IF(IVIS.LE.2) THEN
  703. CALL INDICA(DK1,DK2,IFIS,IPLA,3,BETJEF)
  704. ELSE
  705. CALL INDIC1(DK1,DK2,IFIS,IPLA,3,BETJEF,NECH0,NECH1)
  706. ENDIF
  707. C ********************** CALCUL VISCOPLASTIQUE ***********************
  708. IF (IVIS.EQ.1) THEN
  709. CALL VISPLA(SIGR,SIGF,DSIGT,NSTRS,DEPS1,DEPS2,
  710. & SIGP,SIGRV,DSTRN,D,BETJEF,VISCO)
  711. ENDIF
  712. C ********************************************************************
  713. 100 CONTINUE
  714. IEC1=0
  715. IEC2=0
  716. CRIMAX=0.D0
  717. C
  718. RETURN
  719. END
  720.  
  721.  
  722.  
  723.  
  724.  
  725.  
  726.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales