Télécharger behav1.eso

Retour à la liste

Numérotation des lignes :

  1. C BEHAV1 SOURCE CB215821 16/04/21 21:15:17 8920
  2. SUBROUTINE BEHAV1(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 POST PIC
  9. C Un seul critere de traction: Rankine
  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),SIGRV(4)
  17. DIMENSION V1(4),AC(4),SIGE(4),SIGP(4),SIGR(4),DSTRN(4)
  18. DIMENSION D2FSIG(4,4),DEP(4,4),P1(4,4),D(4,4),DP(4,4),A(4,4)
  19. DIMENSION 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. C
  38. * COMMON /DBETJEF/AA,BETA,COLI,PALF,YOUN,XNU,GFC,GFT,CAR,ETA,TDEF,
  39. * & TCON,DPSTF1,DPSTF2,TETA,PDT,ICT,ICC,IMOD,IVIS,ITR,
  40. * & ISIM,IBB,IGAU,IZON
  41. C
  42. C
  43. IAPEX=0
  44. PRB=1.D-10
  45. PRB2=1.D-6
  46. ITER=1
  47. ITANG=0
  48. IBROY=0
  49. Ft=PALF*COLI
  50. CRIMAX=0.D0
  51. SEQ = 0.D0
  52. SEQ1 = 0.D0
  53. CALL ZERO(SIGE,4,1)
  54. CALL ZERO(V1,4,1)
  55. CALL ZERO(P1,4,4)
  56. CALL ZERO(D,4,4)
  57. C
  58. DO 10 I=1,NSTRS
  59. SIGE(I)=SIGF(I)
  60. 10 CONTINUE
  61. C
  62. C DO 11 I=1,NSTRS
  63. C DO 11 J=1,NSTRS
  64. C D(I,J)=0.D0
  65. C 11 CONTINUE
  66. C
  67. IF (IMOD.EQ.1.OR.IMOD.EQ.3) THEN
  68. AD=YOUN/(1.D0-XNU*XNU)
  69. D(1,1)=AD
  70. D(2,2)=D(1,1)
  71. D(3,3)=AD*(1.D0-XNU)/2.D0
  72. D(1,2)=AD*XNU
  73. D(2,1)=D(1,2)
  74. ENDIF
  75. C
  76. IF (IMOD.EQ.2.OR.IMOD.EQ.4) THEN
  77. ADD=YOUN/((1.D0+XNU)*(1.D0-2.D0*XNU))
  78. D(1,1)=ADD*(1.D0-XNU)
  79. D(2,2)=D(1,1)
  80. D(3,3)=D(1,1)
  81. D(1,2)=ADD*XNU
  82. D(2,1)=D(1,2)
  83. D(1,3)=D(1,2)
  84. D(2,3)=D(1,2)
  85. D(3,1)=D(1,2)
  86. D(3,2)=D(1,2)
  87. D(4,4)=0.5*ADD*(1.D0-2.D0*XNU)
  88. ENDIF
  89. C
  90. C ************ Le point est fissure pour 1ere fois ***********
  91. C
  92. IF (IFIS.EQ.0) THEN
  93. CALL PRINC(SIGF,V1,NSTRS)
  94. DEPS1=0.D0
  95. SIG1=Ft
  96. IF (V1(1).GT.Ft) THEN
  97. IFIS=1
  98. ENDIF
  99. ENDIF
  100. IF (IPLA.EQ.0) THEN
  101. DEPS2=0.D0
  102. SIG2=COLI*AA
  103. ENDIF
  104. C
  105. 8 CONTINUE
  106. C
  107. DO 4 I=1,NSTRS
  108. SIGF(I)=SIGE(I)
  109. 4 CONTINUE
  110. C
  111. CALL PRINC(SIGF,V1,NSTRS)
  112. TETA=V1(4)
  113. DK=DEPS1
  114. DLAM0=0.D0
  115. FCRI0=V1(1)-SIG1
  116. CRIMAX=ABS(100.D0*FCRI0)
  117. C
  118. IF (FCRI0.LT.0.D0) THEN
  119. DO 45 I=1,NSTRS
  120. DO 45 J=1,NSTRS
  121. DEP(I,J)=D(I,J)
  122. 45 CONTINUE
  123. GOTO 100
  124. ENDIF
  125. C
  126. C ************ Traitement du point deja fissure ********************
  127. C
  128. 9 CONTINUE
  129. PI=4.D0*ATAN(1.D0)
  130. PHIC=V1(4)*(PI/180.D0)
  131. COSA=COS(PHIC)
  132. SINA=SIN(PHIC)
  133. C-------------------------------------------------------------------
  134. IF (IMOD.EQ.1.OR.IMOD.EQ.3) THEN
  135. DFSIG(1)=COSA*COSA
  136. DFSIG(2)=SINA*SINA
  137. DFSIG(3)=2.D0*SINA*COSA
  138. DO 12 I=1,NSTRS
  139. DO 12 J=1,NSTRS
  140. P1(I,J)=0.D0
  141. 12 CONTINUE
  142. P1(1,1)=1.D0/2.D0
  143. P1(1,2)=-1.D0/2.D0
  144. P1(2,1)=-1.D0/2.D0
  145. P1(2,2)=1.D0/2.D0
  146. P1(3,3)=2.D0
  147. ENDIF
  148. C-------------------------------------------------------------------
  149. IF (IMOD.EQ.2.OR.IMOD.EQ.4) THEN
  150. DFSIG(1)=COSA*COSA
  151. DFSIG(2)=SINA*SINA
  152. DFSIG(3)=0.D0
  153. DFSIG(4)=2.D0*SINA*COSA
  154. DO 13 I=1,NSTRS
  155. DO 13 J=1,NSTRS
  156. P1(I,J)=0.D0
  157. 13 CONTINUE
  158. P1(1,1)=1.D0/2.D0
  159. P1(1,2)=-1.D0/2.D0
  160. P1(2,1)=-1.D0/2.D0
  161. P1(2,2)=1.D0/2.D0
  162. P1(4,4)=2.D0
  163. ENDIF
  164. C-------------------------------------------------------------------
  165. DO 20 I=1,NSTRS
  166. AC(I)=0.D0
  167. DO 20 J=1,NSTRS
  168. AC(I)=AC(I)+D(I,J)*DFSIG(I)
  169. 20 CONTINUE
  170. C
  171. FDF=0.D0
  172. DO 30 J=1,NSTRS
  173. FDF=FDF+DFSIG(J)*AC(J)
  174. 30 CONTINUE
  175. C--------------- Determination du parametre d'ecrouissage ----------
  176. C
  177. IF(IVIS.LE.2) THEN
  178. CALL UNICOU(DK,PAEC,1,SEQ,BETJEF)
  179. ELSE
  180. CALL UNICO1(DK,PAEC,1,SEQ,BETJEF,NECH0,NECH1)
  181. ENDIF
  182. C
  183. C-------------------------------------------------------------------
  184. DJAC0=-(PAEC+FDF)
  185. C
  186. C ************ Debut des iterations internes ***************
  187. C
  188. 40 CONTINUE
  189. C
  190. C *************** Determination de DK *********************
  191. C
  192. DLAM1=-FCRI0/DJAC0+DLAM0
  193. DK=DEPS1+DLAM1
  194. C IF (DLAM1.LT.0.D0) THEN
  195. C WRITE(*,*)'Dans BEHAV1, DLAMDA est negatif:',DLAM1
  196. C WRITE(*,*)'A l iteration :',ITER
  197. C ENDIF
  198. C
  199. C--------------- Estimation contrainte quivalente ----------
  200. C
  201. IF(IVIS.LE.2) THEN
  202. CALL UNICOU(DK,PAEC,1,SEQ1,BETJEF)
  203. ELSE
  204. CALL UNICO1(DK,PAEC,1,SEQ1,BETJEF,NECH0,NECH1)
  205. ENDIF
  206. C
  207. C ************** Determination de DPHI *********************
  208. C
  209. IF (IMOD.EQ.1.OR.IMOD.EQ.3) THEN
  210. AD1=YOUN/(1.D0-XNU)
  211. TO=SIGF(3)
  212. ENDIF
  213. IF (IMOD.EQ.2.OR.IMOD.EQ.4) THEN
  214. AD1=YOUN/((1.D0+XNU)*(1.D0-2.D0*XNU))
  215. TO=SIGF(4)
  216. ENDIF
  217. C
  218. DPHI=SEQ1-0.5*(SIGE(1)+SIGE(2))+0.5*AD1*DLAM1
  219. C IF (ITER.EQ.1) THEN
  220. C DPHI=0.25*(SIGF(1)-SIGF(2))*(SIGF(1)-SIGF(2))
  221. C DPHI=DPHI+TO*TO
  222. C DPHI=SQRT(DPHI)
  223. C ENDIF
  224. IF (DPHI.LT.0.D0) THEN
  225. C WRITE(*,*)'ATTENTION DPHI NEGATIF'
  226. C WRITE(*,*)'DPHI=',DPHI
  227. ENDIF
  228. C
  229. C--------------- Cas de l'apex -----------------------------
  230. C
  231. IF (ABS(DPHI).LE.10E-10) THEN
  232. IAPEX=1
  233. C WRITE(*,*)'IAPEX ds BEHAV1 =',IAPEX
  234. C WRITE(*,*)'Dans l element',IBB
  235. C WRITE(*,*)'et au point d intégration',IGAU
  236. DO 50 I=1,NSTRS
  237. DO 50 J=1,NSTRS
  238. AI(I,J)=0.D0
  239. 50 CONTINUE
  240. AI(1,1)=0.5
  241. AI(1,2)=0.5
  242. AI(2,1)=AI(1,2)
  243. AI(2,2)=AI(1,1)
  244. IF (IMOD.EQ.1.OR.IMOD.EQ.3) AI(3,3)=0.D0
  245. IF (IMOD.EQ.2.OR.IMOD.EQ.4) AI(3,3)=1.D0
  246. GOTO 75
  247. ENDIF
  248. C
  249. C ************** Mise a jour des contraintes ***************
  250. C
  251. C ---------------- calcul de la matrice A ------------------
  252. C
  253. DO 60 I=1,NSTRS
  254. DO 60 J=1,NSTRS
  255. A(I,J)=0.D0
  256. 60 CONTINUE
  257. C
  258. DG=YOUN/2.D0/(1.D0+XNU)
  259. C
  260. IF (IMOD.EQ.1.OR.IMOD.EQ.3) THEN
  261. A(1,1)=1.D0+(DLAM1*DG)/2.D0/DPHI
  262. A(2,2)=A(1,1)
  263. A(1,2)=-(A(1,1)-1.D0)
  264. A(2,1)=A(1,2)
  265. A(3,3)=1.D0+(DLAM1*DG)/DPHI
  266. ELSE
  267. A(1,1)=1.D0+(DLAM1*DG)/2.D0/DPHI
  268. A(2,2)=A(1,1)
  269. A(1,2)=-(A(1,1)-1.D0)
  270. A(2,1)=A(1,2)
  271. A(3,3)=1.D0
  272. A(4,4)=1.D0+(DLAM1*DG)/DPHI
  273. ENDIF
  274. C
  275. C -------------- invertion de la matrice A -----------------
  276. C
  277. DO 70 I=1,NSTRS
  278. DO 70 J=1,NSTRS
  279. AI(I,J)=A(I,J)
  280. 70 CONTINUE
  281. CALL INVMA2(AI,NSTRS,ISING)
  282. IF (ISING.EQ.1) THEN
  283. WRITE(*,*)'MATRICE AI singuliere ds BEHAV1'
  284. ENDIF
  285. C
  286. C -------------- mise a jour des contraintes ------------
  287. C
  288. 75 CONTINUE
  289. IF (IMOD.EQ.1.OR.IMOD.EQ.3) THEN
  290. VEC1(1)=AD1
  291. VEC1(2)=AD1
  292. VEC1(3)=0.D0
  293. VEC1(4)=0.D0
  294. ENDIF
  295. IF (IMOD.EQ.2.OR.IMOD.EQ.4) THEN
  296. VEC1(1)=AD1
  297. VEC1(2)=AD1
  298. VEC1(3)=AD1*2.D0*XNU
  299. VEC1(4)=0.D0
  300. ENDIF
  301. C
  302. DO 80 I=1,NSTRS
  303. DEPSI(I)=SIGE(I)-0.5*DLAM1*VEC1(I)
  304. 80 CONTINUE
  305. C
  306. DO 90 I=1,NSTRS
  307. SIGF(I)=0.0D+00
  308. DO 90 J=1,NSTRS
  309. SIGF(I)=SIGF(I)+AI(I,J)*DEPSI(J)
  310. 90 CONTINUE
  311. C
  312. C ******** Verification du critere ****************
  313. C
  314. CALL PRINC(SIGF,V1,NSTRS)
  315. FCRI1=V1(1)-SEQ1
  316. C
  317. IF (IBROY.EQ.0.AND.ABS(FCRI1).GE.CRIMAX) THEN
  318. C WRITE(*,*)'****************************************'
  319. C WRITE(*,*)'LE RESIDU DIVERGE AVEC BROYDEN'
  320. C WRITE(*,*)'on passe donc a la secante'
  321. C WRITE(*,*)'Dans l element',IBB
  322. C WRITE(*,*)'et au point d intégration',IGAU
  323. C WRITE(*,*)'CRIMAX=',CRIMAX
  324. C WRITE(*,*)'****************************************'
  325. ITER=ITR
  326. ENDIF
  327. C
  328. C ******* Compteur sur la methode de resolution ****
  329. C
  330. IF (IBROY.EQ.0.AND.ITER.EQ.ITR) THEN
  331. IBROY=1
  332. ITANG=1
  333. ITER=1
  334. IAPEX=0
  335. GOTO 8
  336. ENDIF
  337. C
  338. C ******* non convergence **************************
  339. C
  340. IF (ABS(FCRI1).GT.PRB.AND.ITER.LT.ITR) THEN
  341. IF (IBROY.EQ.0) THEN
  342. DJAC1=(FCRI0-FCRI1)/(DLAM0-DLAM1)
  343. DLAM0=DLAM1
  344. DJAC0=DJAC1
  345. FCRI0=FCRI1
  346. ITER=ITER+1
  347. IF (ITER.GE.(ITR-1)) THEN
  348. C WRITE(*,*)'***********************'
  349. C WRITE(*,*)'BROYDEN n a pas aboutit'
  350. C WRITE(*,*)'Dans l element',IBB
  351. C WRITE(*,*)'et au point d intégration',IGAU
  352. C WRITE(*,*)'ITER=',ITER
  353. C WRITE(*,*)'FCRI=',FCRI1
  354. C WRITE(*,*)'***********************'
  355. ENDIF
  356. GOTO 40
  357. ENDIF
  358. IF (IBROY.EQ.1.AND.ITANG.EQ.1) THEN
  359. DLAM0=DLAM1
  360. FCRI0=FCRI1
  361. CALL PRINC(SIGF,V1,NSTRS)
  362. ITER=ITER+1
  363. IF (ITER.GE.(ITR-5)) THEN
  364. WRITE(*,*)'ITER=',ITER
  365. WRITE(*,*)'FCRI=',FCRI1
  366. ENDIF
  367. GOTO 9
  368. ENDIF
  369. ENDIF
  370. IF (ITER.GE.ITR.AND.ABS(FCRI1).GT.PRB2) THEN
  371. WRITE(*,*)'NON CONVERGENCE INTERNE dans BEHAV1'
  372. WRITE(*,*)'FCRI=',FCRI1
  373. WRITE(*,*)'Dans l element',IBB
  374. WRITE(*,*)'et au point d intégration',IGAU
  375. WRITE(*,*)'IPLA=',IPLA
  376. WRITE(*,*)'IFIS=',IFIS
  377. C STOP
  378. ENDIF
  379. C
  380. C ************* Calcul de la DEP ****************
  381. C
  382. IF (IAPEX.EQ.0) THEN
  383. CALL DERI2(D2FSIG,DPHI,P1,SIGF,NSTRS,BETJEF)
  384. CALL LADEP(SIGF,DEP,PAEC,DPHI,NSTRS,IFOU,
  385. & DFSIG,D2FSIG,DLAM1,D,DP,BETJEF)
  386. ELSE
  387. C WRITE(*,*)'ELAS',IAPEX
  388. DO 95 I=1,NSTRS
  389. DO 95 J=1,NSTRS
  390. AH(I,J)=0.D0
  391. DO 95 K=1,NSTRS
  392. AH(I,J)=AH(I,J)+AI(I,K)*D(K,J)
  393. 95 CONTINUE
  394. DFSIG(1)=0.5*SQRT(2.D0)
  395. DFSIG(2)=0.5*SQRT(2.D0)
  396. DFSIG(3)=0.D0
  397. DFSIG(4)=0.D0
  398. CALL CREDEP(AH,DFSIG,PAEC,NSTRS,DEP,BETJEF)
  399. ENDIF
  400. IAPEX=0
  401. CRIMAX=0.D0
  402. C
  403. C ***********************************************
  404. C
  405. IF (DK.LT.0.D0) THEN
  406. WRITE(*,*)'Dans BEHAV1, DLAMDA est negatif:',DLAM1
  407. WRITE(*,*)'A l iteration :',ITER
  408. C STOP
  409. ENDIF
  410. DEPS1=DK
  411. SIG1=V1(1)
  412. IF (V1(1).LT.0.D0) V1(1)=0.D0
  413. CALL INDICA(DK,0.0D0,IFIS,IPLA,1,BETJEF)
  414. IF(IVIS.LE.2) THEN
  415. CALL INDICA(DK,0.0D0,IFIS,IPLA,1,BETJEF)
  416. ELSE
  417. CALL INDIC1(DK,0.0D0,IFIS,IPLA,1,BETJEF,NECH0,NECH1)
  418. ENDIF
  419. C ********************** CALCUL VISCOPLASTIQUE ***********************
  420. IF (IVIS.EQ.1) THEN
  421. CALL VISPLA(SIGR,SIGF,DSIGT,NSTRS,DEPS1,DEPS2,
  422. & SIGP,SIGRV,DSTRN,D,BETJEF,VISCO)
  423. ENDIF
  424. C ********************************************************************
  425. 100 CONTINUE
  426. C
  427. RETURN
  428. END
  429.  
  430.  
  431.  
  432.  
  433.  
  434.  
  435.  
  436.  
  437.  

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