Télécharger behav2.eso

Retour à la liste

Numérotation des lignes :

  1. C BEHAV2 SOURCE CB215821 16/04/21 21:15:18 8920
  2. SUBROUTINE BEHAV2(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 Un seul critere de compression: Von Mises ou Drucker-Prager
  10. C
  11. C ==================================================================
  12. C CE SOUS-PROGRAMME EST APPELE DANS "BONE".
  13. C
  14. IMPLICIT INTEGER(I-N)
  15. IMPLICIT REAL*8(A-H,O-Z)
  16. DIMENSION DFSIG(4),SIGF(4),DSIGT(4),VEC1(4),DEPSI(4),DSTRN(4)
  17. DIMENSION SIGR(4),AC(4),SIGE(4),SIGRV(4),SIGP(4)
  18. DIMENSION D(4,4),DEP(4,4),D2FSIG(4,4),P2(4,4),DP(4,4)
  19. DIMENSION A(4,4),AI(4,4),AH(4,4)
  20. *
  21. SEGMENT BETJEF
  22. REAL*8 AA,BETA,COLI,PALF,YOUN,XNU,GFC,GFT,CAR,ETA,TDEF,
  23. & TCON,DPSTF1,DPSTF2,TETA,PDT,TP0
  24. INTEGER ICT,ICC,IMOD,IVIS,ITR,
  25. & ISIM,IBB,IGAU,IZON
  26. ENDSEGMENT
  27. *
  28. SEGMENT VISCO
  29. REAL*8 DPSTV1,DPSTV2,SIGV1,SIGV2,ENDV
  30. ENDSEGMENT
  31. SEGMENT NECH0
  32. REAL*8 DT,DC,ALFG,S0,ENDO
  33. ENDSEGMENT
  34. SEGMENT NECH1
  35. REAL*8 ENDL
  36. ENDSEGMENT
  37. *
  38. *
  39. C
  40. * COMMON /DBETJEF/AA,BETA,COLI,PALF,YOUN,XNU,GFC,GFT,CAR,ETA,TDEF,
  41. * & TCON,DPSTF1,DPSTF2,TETA,PDT,ICT,ICC,IMOD,IVIS,ITR,
  42. * & ISIM,IBB,IGAU,IZON
  43. C
  44. IAPEX=0
  45. PRB=1.D-10
  46. PRB2=1.D-6
  47. ITER=1
  48. IBROY=0
  49. ITANG=0
  50. Rb=COLI
  51. Ft=PALF*COLI
  52. IEC2=0
  53. CRIMAX=0.D0
  54. SEQ = 0.D0
  55. SEQ2 = 0.D0
  56. CALL ZERO(SIGE,4,1)
  57. ***** CALL ZERO(V1,4,1)
  58. CALL ZERO(P2,4,4)
  59. CALL ZERO(D,4,4)
  60. C
  61. DO 10 I=1,NSTRS
  62. SIGE(I)=SIGF(I)
  63. 10 CONTINUE
  64. C
  65. IF (IMOD.EQ.1.OR.IMOD.EQ.3) THEN
  66. AD=YOUN/(1.D0-XNU*XNU)
  67. D(1,1)=AD
  68. D(2,2)=D(1,1)
  69. D(3,3)=AD*(1.D0-XNU)/2.D0
  70. D(1,2)=AD*XNU
  71. D(2,1)=D(1,2)
  72. ENDIF
  73. C
  74. IF (IMOD.EQ.2.OR.IMOD.EQ.4) THEN
  75. ADD=YOUN/((1.D0+XNU)*(1.D0-2.D0*XNU))
  76. D(1,1)=ADD*(1.D0-XNU)
  77. D(2,2)=D(1,1)
  78. D(3,3)=D(1,1)
  79. D(1,2)=ADD*XNU
  80. D(2,1)=D(1,2)
  81. D(1,3)=D(1,2)
  82. D(2,3)=D(1,2)
  83. D(3,1)=D(1,2)
  84. D(3,2)=D(1,2)
  85. D(4,4)=0.5*ADD*(1.D0-2.D0*XNU)
  86. ENDIF
  87. C
  88. C ************ Le point est comprime pour 1ere fois ***********
  89. C
  90. IF (IPLA.EQ.0) THEN
  91. IF(IMOD.EQ.1.OR.IMOD.EQ.2) THEN
  92. CALL CRIVON(SIGF,SEQ,NSTRS,BETJEF)
  93. ENDIF
  94. IF(IMOD.EQ.3.OR.IMOD.EQ.4) THEN
  95. CALL CRIDRU(SIGF,SEQ,NSTRS,BETJEF)
  96. ENDIF
  97. DEPS2=0.D0
  98. SIG2=Rb*AA
  99. IF (SEQ.GT.SIG2) THEN
  100. IPLA=1
  101. ENDIF
  102. ENDIF
  103. IF (IFIS.EQ.0) THEN
  104. DEPS1=0.D0
  105. SIG1=Ft
  106. ENDIF
  107. C
  108. 8 CONTINUE
  109. C
  110. IF(IMOD.EQ.1.OR.IMOD.EQ.2) THEN
  111. CALL CRIVON(SIGE,SEQ,NSTRS,BETJEF)
  112. ENDIF
  113. IF(IMOD.EQ.3.OR.IMOD.EQ.4) THEN
  114. CALL CRIDRU(SIGE,SEQ,NSTRS,BETJEF)
  115. ENDIF
  116. FCRIC0=SEQ-SIG2
  117. CRIMAX=ABS(100.D0*FCRIC0)
  118. IF (FCRIC0.LE.0.D0) THEN
  119. IEC2=1
  120. ENDIF
  121. DK=DEPS2
  122. DLAM0=0.D0
  123. C
  124. 12 CONTINUE
  125. C
  126. DO 15 I=1,NSTRS
  127. SIGF(I)=SIGE(I)
  128. 15 CONTINUE
  129. C
  130. IF (IEC2.EQ.1) THEN
  131. DO 16 I=1,NSTRS
  132. DO 16 J=1,NSTRS
  133. DEP(I,J)=D(I,J)
  134. 16 CONTINUE
  135. GOTO 100
  136. ENDIF
  137. C
  138. 9 CONTINUE
  139. C
  140. IF(IVIS.LE.2) THEN
  141. CALL UNICOU(DK,PAEC,2,SEQ2,BETJEF)
  142. ELSE
  143. CALL UNICO1(DK,PAEC,2,SEQ2,BETJEF,NECH0,NECH1)
  144. ENDIF
  145. C
  146. C ************ Preparation des criteres ********************
  147. C
  148. C ------------------ Von Mises -----------------------------
  149. IF (IMOD.EQ.1.OR.IMOD.EQ.2) THEN
  150. A1=0.D0
  151. A2=1.D0
  152. ENDIF
  153. C ----------------- Drucker Prager -------------------------
  154. IF (IMOD.EQ.3.OR.IMOD.EQ.4) THEN
  155. A1=(BETA-1.D0)/(2*BETA-1.D0)
  156. A2=BETA/(2*BETA-1.D0)
  157. ENDIF
  158. C ----------------------------------------------------------
  159. IF (IMOD.EQ.1.OR.IMOD.EQ.3) THEN
  160. SX=SIGF(1)
  161. SY=SIGF(2)
  162. SXY=SIGF(3)
  163. R1=(SX-SY)*(SX-SY)+SX*SX+SY*SY+6*SXY*SXY
  164. R=SQRT(0.5*R1)
  165. IF (ABS(R).LE.1.D-10) THEN
  166. WRITE(*,*)'INDETERMINATION de DFSIG ds BEHAV2'
  167. STOP
  168. ENDIF
  169. DFSIG(1)=(2.D0*SX-SY)/2.D0/R/A2+A1/A2
  170. DFSIG(2)=(2.D0*SY-SX)/2.D0/R/A2+A1/A2
  171. DFSIG(3)=3.D0*SXY/R/A2
  172. P2(1,1)=2.D0/A2/A2
  173. P2(1,2)=-1.D0/A2/A2
  174. P2(1,3)=0.D0/A2/A2
  175. P2(2,1)=-1.D0/A2/A2
  176. P2(2,2)=2.D0/A2/A2
  177. P2(2,3)=0.D0/A2/A2
  178. P2(3,1)=0.D0/A2/A2
  179. P2(3,2)=0.D0/A2/A2
  180. P2(3,3)=6.D0/A2/A2
  181. ENDIF
  182. C
  183. IF (IMOD.EQ.2.OR.IMOD.EQ.4) THEN
  184. SX=SIGF(1)
  185. SY=SIGF(2)
  186. SZ=SIGF(3)
  187. SXY=SIGF(4)
  188. TX=(SX-SY)*(SX-SY)
  189. TY=(SX-SZ)*(SX-SZ)
  190. TZ=(SY-SZ)*(SY-SZ)
  191. R1=TX+TY+TZ+6*SXY*SXY
  192. R=SQRT(0.5*R1)
  193. IF (ABS(R).LE.1.D-10) THEN
  194. WRITE(*,*)'INDETERMINATION de DFSIG ds BEHAV2'
  195. ENDIF
  196. DFSIG(1)=(2.D0*SX-SY-SZ)/2.D0/R/A2+A1/A2
  197. DFSIG(2)=(2.D0*SY-SX-SZ)/2.D0/R/A2+A1/A2
  198. DFSIG(3)=(2.D0*SZ-SX-SY)/2.D0/R/A2+A1/A2
  199. DFSIG(4)=3.D0*SXY/R/A2
  200. P2(1,1)=2.D0/A2/A2
  201. P2(1,2)=-1.D0/A2/A2
  202. P2(1,3)=-1.D0/A2/A2
  203. P2(1,4)=0.D0/A2/A2
  204. P2(2,1)=-1.D0/A2/A2
  205. P2(2,2)=2.D0/A2/A2
  206. P2(2,3)=-1.D0/A2/A2
  207. P2(2,4)=0.D0/A2/A2
  208. P2(3,1)=-1.D0/A2/A2
  209. P2(3,2)=-1.D0/A2/A2
  210. P2(3,3)=2.D0/A2/A2
  211. P2(3,4)=0.D0/A2/A2
  212. P2(4,1)=0.D0/A2/A2
  213. P2(4,2)=0.D0/A2/A2
  214. P2(4,3)=0.D0/A2/A2
  215. P2(4,4)=6.D0/A2/A2
  216. ENDIF
  217. C
  218. DO 25 I=1,NSTRS
  219. AC(I)=0.D0
  220. DO 25 J=1,NSTRS
  221. AC(I)=AC(I)+D(I,J)*DFSIG(J)
  222. 25 CONTINUE
  223. C
  224. C ************ Preparation du jacobien ***********************
  225. C
  226. FDF=0.D0
  227. DO 30 J=1,NSTRS
  228. FDF=FDF+DFSIG(J)*AC(J)
  229. 30 CONTINUE
  230. C
  231. DJAC0=-(FDF+PAEC)
  232. C
  233. C ************ Debut des iterations internes *******************
  234. C
  235. 40 CONTINUE
  236. C
  237. C *************** Determination de DK et DLAM ******************
  238. C
  239. C
  240. DLAM1=DLAM0-FCRIC0/DJAC0
  241. C
  242. IF (DLAM1.LT.0.D0) IEC2=1
  243. IF (IEC2.EQ.1) THEN
  244. C WRITE(*,*)'Dans BEHAV2, DLAMDA est negatif:',DLAM1
  245. C WRITE(*,*)'A l iteration :',ITER
  246. GOTO 12
  247. ENDIF
  248. C
  249. DK=DEPS2+DLAM1
  250. C
  251. IF(IVIS.LE.2) THEN
  252. CALL UNICOU(DK,PAEC,2,SEQ2,BETJEF)
  253. ELSE
  254. CALL UNICO1(DK,PAEC,2,SEQ2,BETJEF,NECH0,NECH1)
  255. ENDIF
  256. C
  257. C ************** Determination de DPHI2 *********************
  258. C
  259. IF (IMOD.EQ.1.OR.IMOD.EQ.2) THEN
  260. DPHI2=SEQ2
  261. ENDIF
  262. IF (IMOD.EQ.3) THEN
  263. SOM=SIGF(1)+SIGF(2)+SIGF(3)
  264. DPHI2=SEQ2-A1*SOM/A2
  265. ENDIF
  266. IF (IMOD.EQ.4) THEN
  267. SOM=SIGE(1)+SIGE(2)+SIGE(3)
  268. DPHI2=A1*SOM/A2-A1*A1*DLAM1*3.D0*ADD*(1.D0+XNU)/A2/A2
  269. DPHI2=SEQ2-DPHI2
  270. ENDIF
  271. C IF (DPHI2.LT.0.D0) THEN
  272. C WRITE(*,*)'ATTENTION DPHI2 ds BEHAV2 NEGATIF'
  273. C WRITE(*,*)'DPHI2=',DPHI2
  274. C ENDIF
  275. C
  276. C---------- Cas du critere reduit a un point ----------------
  277. C
  278. IF (ABS(DPHI2).LE.10E-10) THEN
  279. IAPEX=1
  280. C WRITE(*,*)'IAPEX ds BEHAV2 =',IAPEX
  281. DO 50 I=1,NSTRS
  282. DO 50 J=1,NSTRS
  283. AI(I,J)=0.D0
  284. 50 CONTINUE
  285. IF (IMOD.EQ.2.OR.IMOD.EQ.4) THEN
  286. DO 53 I=1,3
  287. DO 53 J=1,3
  288. AI(I,J)=1.D0/3.D0
  289. 53 CONTINUE
  290. ENDIF
  291. GOTO 75
  292. ENDIF
  293. C
  294. C ************** Mise a jour des contraintes ***************
  295. C
  296. C ---------------- calcul de la matrice A ------------------
  297. C
  298. DO 60 I=1,NSTRS
  299. DO 60 J=1,NSTRS
  300. A(I,J)=0.D0
  301. 60 CONTINUE
  302. C
  303. DG=YOUN/2.D0/(1.D0+XNU)
  304. C
  305. IF (IMOD.EQ.1.OR.IMOD.EQ.3) THEN
  306. DW1=AD*(2.D0-XNU)
  307. DW2=AD*(2.D0*XNU-1.D0)
  308. A(1,1)=1.D0+(DLAM1*DW1)/2.D0/DPHI2/A2/A2
  309. A(2,2)=A(1,1)
  310. A(1,2)=(DLAM1*DW2)/2.D0/DPHI2/A2/A2
  311. A(2,1)=A(1,2)
  312. A(3,3)=1.D0+(DLAM1*3.D0*DG)/DPHI2/A2/A2
  313. ENDIF
  314. IF (IMOD.EQ.2.OR.IMOD.EQ.4) THEN
  315. A(1,1)=1.D0+(DLAM1*2.D0*DG)/DPHI2/A2/A2
  316. A(2,2)=A(1,1)
  317. A(3,3)=1.D0+(DLAM1*2.D0*DG)/DPHI2/A2/A2
  318. A(1,2)=-(DLAM1*DG)/DPHI2/A2/A2
  319. A(2,1)=A(1,2)
  320. A(1,3)=-(DLAM1*DG)/DPHI2/A2/A2
  321. A(3,1)=A(1,3)
  322. A(3,2)=A(1,3)
  323. A(2,3)=A(1,3)
  324. A(4,4)=1.D0+(DLAM1*DG*3.D0)/DPHI2/A2/A2
  325. ENDIF
  326. C
  327. C -------------- invertion de la matrice A -----------------
  328. C
  329. DO 70 I=1,NSTRS
  330. DO 70 J=1,NSTRS
  331. AI(I,J)=A(I,J)
  332. 70 CONTINUE
  333. CALL INVMA2(AI,NSTRS,ISING)
  334. IF (ISING.EQ.1) THEN
  335. WRITE(*,*)'MATRICE AI singuliere ds BEHAV2'
  336. WRITE(*,*)'NSTRS=',NSTRS
  337. ENDIF
  338. C
  339. C -------------- mise a jour des contraintes --------------
  340. C
  341. 75 CONTINUE
  342. IF(IMOD.EQ.3) THEN
  343. VEC1(1)=AD*(1.D0+XNU)
  344. VEC1(2)=AD*(1.D0+XNU)
  345. VEC1(3)=0.D0
  346. VEC1(4)=0.D0
  347. ENDIF
  348. IF(IMOD.EQ.4) THEN
  349. VEC1(1)=ADD*(1.D0+XNU)
  350. VEC1(2)=ADD*(1.D0+XNU)
  351. VEC1(3)=ADD*(1.D0+XNU)*2.D0*XNU
  352. VEC1(4)=0.D0
  353. ENDIF
  354. IF(IMOD.EQ.1.OR.IMOD.EQ.2) THEN
  355. VEC1(1)=0.D0
  356. VEC1(2)=0.D0
  357. VEC1(3)=0.D0
  358. VEC1(4)=0.D0
  359. ENDIF
  360. C
  361. DO 80 I=1,NSTRS
  362. DEPSI(I)=SIGE(I)-A1*DLAM1*VEC1(I)/A2
  363. 80 CONTINUE
  364. C
  365. DO 90 I=1,NSTRS
  366. SIGF(I)=0.0D+00
  367. DO 90 J=1,NSTRS
  368. SIGF(I)=SIGF(I)+AI(I,J)*DEPSI(J)
  369. 90 CONTINUE
  370. C
  371. C ************ Verification du critere ********************
  372. C
  373. IF(IMOD.EQ.1.OR.IMOD.EQ.2) THEN
  374. CALL CRIVON(SIGF,SEQ,NSTRS,BETJEF)
  375. ENDIF
  376. IF(IMOD.EQ.3.OR.IMOD.EQ.4) THEN
  377. CALL CRIDRU(SIGF,SEQ,NSTRS,BETJEF)
  378. ENDIF
  379. FCRIC1=SEQ-SEQ2
  380. C
  381. IF (IBROY.EQ.0.AND.ABS(FCRIC1).GE.CRIMAX) THEN
  382. C WRITE(*,*)'****************************************'
  383. C WRITE(*,*)'LE RESIDU DIVERGE AVEC BROYDEN'
  384. C WRITE(*,*)'on passe donc a la secante'
  385. C WRITE(*,*)'Dans l element',IBB
  386. C WRITE(*,*)'et au point d intégration',IGAU
  387. C WRITE(*,*)'CRIMAX=',CRIMAX
  388. C WRITE(*,*)'****************************************'
  389. ITER=ITR
  390. ENDIF
  391. C
  392. C ******* Compteur sur la methode de resolution ****
  393. C
  394. IF (IBROY.EQ.0.AND.ITER.EQ.ITR) THEN
  395. IBROY=1
  396. ITANG=1
  397. IAPEX=0
  398. IEC2=0
  399. ITER=1
  400. GOTO 8
  401. ENDIF
  402. C
  403. C ******* non convergence **************************
  404. C
  405. IF (ABS(FCRIC1).GT.PRB.AND.ITER.LT.ITR) THEN
  406. IF (IBROY.EQ.0) THEN
  407. DJAC0=(FCRIC0-FCRIC1)/(DLAM0-DLAM1)
  408. DLAM0=DLAM1
  409. FCRIC0=FCRIC1
  410. ITER=ITER+1
  411. IF (ITER.GE.(ITR-1)) THEN
  412. C WRITE(*,*)'***********************'
  413. C WRITE(*,*)'BROYDEN n a pas aboutit'
  414. C WRITE(*,*)'Dans l element',IBB
  415. C WRITE(*,*)'et au point d intégration',IGAU
  416. C WRITE(*,*)'ITER=',ITER
  417. C WRITE(*,*)'FCRI=',FCRIC1
  418. C WRITE(*,*)'***********************'
  419. ENDIF
  420. GOTO 40
  421. ENDIF
  422. IF (IBROY.EQ.1.AND.ITANG.EQ.1) THEN
  423. DLAM0=DLAM1
  424. FCRIC0=FCRIC1
  425. ITER=ITER+1
  426. IF (ITER.GE.(ITR-5)) THEN
  427. C WRITE(*,*)'ITER=',ITER
  428. C WRITE(*,*)'FCRI=',FCRIC1
  429. ENDIF
  430. GOTO 9
  431. ENDIF
  432. ENDIF
  433. IF (ITER.GE.ITR.AND.ABS(FCRIC1).GT.PRB2) THEN
  434. WRITE(*,*)'NON CONVERGENCE INTERNE dans BEHAV2'
  435. WRITE(*,*)'FCRI=',FCRIC1
  436. WRITE(*,*)'Dans l element',IBB
  437. WRITE(*,*)'et au point d intégration',IGAU
  438. WRITE(*,*)'IPLA=',IPLA
  439. WRITE(*,*)'IFIS=',IFIS
  440. C STOP
  441. ENDIF
  442. C
  443. C ************* Calcul de la DEP ****************
  444. C
  445. IF (IAPEX.EQ.0) THEN
  446. CALL DERI2(D2FSIG,DPHI2,P2,SIGF,NSTRS,BETJEF)
  447. CALL LADEP(SIGF,DEP,PAEC,DPHI,NSTRS,IFOU,
  448. & DFSIG,D2FSIG,DLAM1,D,DP,BETJEF)
  449. ELSE IF (IAPEX.EQ.1) THEN
  450. C WRITE(*,*)'ELAS2',IAPEX
  451. DO 96 I=1,NSTRS
  452. DO 96 J=1,NSTRS
  453. DEP(I,J)=0.D0
  454. DO 96 K=1,NSTRS
  455. DEP(I,J)=DEP(I,J)+AI(I,K)*D(K,J)
  456. 96 CONTINUE
  457. ENDIF
  458. IAPEX=0
  459. CRIMAX=0.D0
  460. * DO 97 I=1,NSTRS
  461. * DO 97 J=1,NSTRS
  462. * DEP(I,J)=D(I,J)
  463. * 97 CONTINUE
  464. C
  465. C ***********************************************
  466. C
  467. DEPS2=DK
  468. SIG2=SEQ2
  469. IF(IVIS.LE.2) THEN
  470. CALL INDICA(0.D0,DK,IFIS,IPLA,2,BETJEF)
  471. ELSE
  472. CALL INDIC1(0.D0,DK,IFIS,IPLA,2,BETJEF,NECH0,NECH1)
  473. ENDIF
  474. C ********************** CALCUL VISCOPLASTIQUE ***********************
  475. IF (IVIS.EQ.1) THEN
  476. CALL VISPLA(SIGR,SIGF,DSIGT,NSTRS,DEPS1,DEPS2,
  477. & SIGP,SIGRV,DSTRN,D,BETJEF,VISCO)
  478. ENDIF
  479. C ********************************************************************
  480. 100 CONTINUE
  481. IEC2=0
  482. C
  483. RETURN
  484. END
  485.  
  486.  
  487.  
  488.  
  489.  
  490.  
  491.  
  492.  
  493.  

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