Télécharger acti3.eso

Retour à la liste

Numérotation des lignes :

  1. C ACTI3 SOURCE CB215821 16/04/21 21:15:04 8920
  2. C ACTI3
  3. SUBROUTINE ACTI3(SIG0,SIGF,D,NSTRSS,BETINSA)
  4. C
  5. C ==================================================================
  6. C Deux criteres de traction compression: DRUCKER PRAGER
  7. C 3D
  8. C ==================================================================
  9. C CE SOUS-PROGRAMME EST APPELE DANS "BONE3D".
  10. C
  11. IMPLICIT INTEGER(I-N)
  12. IMPLICIT REAL*8(A-H,O-Z)
  13. DIMENSION SIG0(NSTRSS),SIGF(NSTRSS)
  14. DIMENSION DFSIG1(6),DFSIG2(6),VEC1(6),VEC2(6),SIGE(6)
  15. DIMENSION AC1(6),AC2(6)
  16. DIMENSION DJAC0(2,2),DJI(4,4),A(6,6),AI(6,6),AIM(6,6),D(6,6)
  17. DIMENSION FCRI0(2),FCRI1(2),DEPSI(6),DLAM1(2),DLAM0(2)
  18. C
  19. SEGMENT BETINSA
  20. REAL*8 RT,RC,YOUN,XNU,GFT,GFC,CAR
  21. REAL*8 DKT,DKC,SEQT,SEQC,ENDT,ENDC
  22. INTEGER IFIS,IPLA,IBB,IGAU
  23. ENDSEGMENT
  24. C
  25. IAPEX=0
  26. PRB=1.D-5
  27. PRB2=1.D-2
  28. ITER=1
  29. ITR = 1500
  30. IET=0
  31. IEC=0
  32. IBROY=0
  33. ITANG=0
  34. CRIMAX=0.D0
  35. SEQ = 0.D0
  36. SEQ1 = 0.D0
  37. SEQ2 = 0.D0
  38. CALL ZERO(SIGE,6,1)
  39. C
  40. DO 10 I=1,NSTRSS
  41. SIGE(I)=SIGF(I)
  42. 10 CONTINUE
  43. C
  44. C ************ Le point est fissure pour 1ere fois ***********
  45. C
  46. 12 CONTINUE
  47. CALL DRUTRA(SIGF,SEQTE,BETINSA)
  48. CALL DRUCOM(SIGF,SEQCE,BETINSA)
  49. FCRI0(1) = SEQTE - SEQT
  50. FCRI0(2) = SEQCE - SEQC
  51. CRIMAX=ABS(100.D0*FCRI0(1))
  52. IF (FCRI0(1).LT.0.D0.AND.FCRI0(2).LT.0.D0) THEN
  53. IET=1
  54. IEC=1
  55. ENDIF
  56. IF (FCRI0(1).LT.0.D0.AND.FCRI0(2).GE.0.D0) THEN
  57. IET=1
  58. IEC=0
  59. ENDIF
  60. IF (FCRI0(1).GE.0.D0.AND.FCRI0(2).LT.0.D0) THEN
  61. IET=0
  62. IEC=1
  63. ENDIF
  64. C
  65. 15 CONTINUE
  66. C
  67. DO 16 I=1,NSTRSS
  68. SIGF(I)=SIGE(I)
  69. 16 CONTINUE
  70. C
  71. IF (IET.EQ.1.AND.IEC.EQ.1) THEN
  72. GOTO 100
  73. ENDIF
  74. C
  75. IF (IET.EQ.1.AND.IEC.EQ.0) THEN
  76. CALL ACTI2(SIG0,SIGF,D,NSTRSS,BETINSA)
  77. GOTO 100
  78. ENDIF
  79. C
  80. IF (IET.EQ.0.AND.IEC.EQ.1) THEN
  81. CALL ACTI1(SIG0,SIGF,D,NSTRSS,BETINSA)
  82. GOTO 100
  83. ENDIF
  84. C
  85. DK1=DKT
  86. DK2=DKC
  87. DLAM0(1)=0.D0
  88. DLAM0(2)=0.D0
  89. C
  90. 18 CONTINUE
  91. C
  92. C ************ CALCUL DU JACOBIEN INITIAL ********************
  93. C
  94. C ---------------- direction de traction ------------------
  95. C
  96. CALL DRUTR1(SIGF,DFSIG1,BETINSA)
  97. C
  98. DO 20 I=1,NSTRSS
  99. AC1(I)=0.D0
  100. DO 20 J=1,NSTRSS
  101. AC1(I)=AC1(I)+D(I,J)*DFSIG1(J)
  102. 20 CONTINUE
  103. C
  104. C ---------------- direction de compression ----------------
  105. C
  106. CALL DRUCO1(SIGF,DFSIG2,BETINSA)
  107. C
  108. DO 25 I=1,NSTRSS
  109. AC2(I)=0.D0
  110. DO 25 J=1,NSTRSS
  111. AC2(I)=AC2(I)+D(I,J)*DFSIG2(J)
  112. 25 CONTINUE
  113. C------------------------------------------------------------
  114. F1DF1=0.D0
  115. DO 30 J=1,NSTRSS
  116. F1DF1=F1DF1+DFSIG1(J)*AC1(J)
  117. 30 CONTINUE
  118. F1DF2=0.D0
  119. DO 31 J=1,NSTRSS
  120. F1DF2=F1DF2+DFSIG1(J)*AC2(J)
  121. 31 CONTINUE
  122. F2DF2=0.D0
  123. DO 32 J=1,NSTRSS
  124. F2DF2=F2DF2+DFSIG2(J)*AC2(J)
  125. 32 CONTINUE
  126. F2DF1=0.D0
  127. DO 33 J=1,NSTRSS
  128. F2DF1=F2DF1+DFSIG2(J)*AC1(J)
  129. 33 CONTINUE
  130. C
  131. CALL ENDAME(1,BETINSA)
  132. CALL FORECR(DK1,PAECT,1,SEQ,BETINSA)
  133. CALL ENDAME(2,BETINSA)
  134. CALL FORECR(DK2,PAECC,2,SEQ,BETINSA)
  135. C
  136. DJAC0(1,1)=-(F1DF1+PAECT)
  137. DJAC0(1,2)=-(F1DF2)
  138. DJAC0(2,2)=-(F2DF2+PAECC)
  139. DJAC0(2,1)=-(F2DF1)
  140. C
  141. DO 43 I=1,2
  142. DO 43 J=1,2
  143. DJI(I,J)=DJAC0(I,J)
  144. 43 CONTINUE
  145. C
  146. CALL INVMA2(DJI,2,ISING)
  147. IF (ISING.EQ.1) THEN
  148. WRITE(*,*)'MATRICE DJI singuliere ds ACTI3'
  149. ENDIF
  150. C
  151. C ************ DEBUT ITERATION INTERNES *******************
  152. C
  153. 40 CONTINUE
  154. C
  155. C *************** Determination de DK et DLAM ******************
  156. C
  157. C
  158. C
  159. DLAM1(1)=DLAM0(1)-DJI(1,1)*FCRI0(1)-DJI(1,2)*FCRI0(2)
  160. DLAM1(2)=DLAM0(2)-DJI(2,1)*FCRI0(1)-DJI(2,2)*FCRI0(2)
  161. C
  162. IF (DLAM1(1).LE.0.D0) IET=1
  163. IF (DLAM1(2).LE.0.D0) IEC=1
  164. IF (IET.EQ.1.OR.IEC.EQ.1) THEN
  165. C WRITE(*,*)'Dans ACTI3, DLAMDA1 est negatif:',DLAM1(1)
  166. C WRITE(*,*)'Dans ACTI3, DLAMDA2 est negatif:',DLAM1(2)
  167. C WRITE(*,*)'A l iteration :',ITER
  168. GOTO 15
  169. ENDIF
  170. C
  171. DK1=DKT+DLAM1(1)
  172. DK2=DKC+DLAM1(2)
  173. C
  174. CALL ENDAME(1,BETINSA)
  175. CALL FORECR(DK1,PAECT,1,SEQ1,BETINSA)
  176. CALL ENDAME(2,BETINSA)
  177. CALL FORECR(DK2,PAECC,2,SEQ2,BETINSA)
  178. C
  179. C ************** Determination de DPHI1 et 2 ******************
  180. C
  181. CALL DRUTR2(SIGE,SEQ1,DPHI1,DLAM1(1),VEC1,BETINSA)
  182. CALL DRUCO2(SIGE,SEQ2,DPHI2,DLAM1(2),VEC2,BETINSA)
  183. C
  184. C--------------- Cas de l'apex -----------------------------
  185. C
  186. IF (ABS(DPHI1).LE.10E-10.AND.
  187. *ABS(DPHI2).GT.10E-10) THEN
  188. IAPEX=1
  189. C WRITE(*,*)'IAPEX ds ACTI3 =',IAPEX
  190. C WRITE(*,*)'Dans l element',IBB
  191. C WRITE(*,*)'et au point d intégration',IGAU
  192. DO 50 I=1,NSTRSS
  193. DO 50 J=1,NSTRSS
  194. AI(I,J)=0.D0
  195. 50 CONTINUE
  196. AI(1,1)=1./3.
  197. AI(1,2)=AI(1,1)
  198. AI(1,2)=AI(1,1)
  199. AI(1,3)=AI(1,1)
  200. AI(2,1)=AI(1,1)
  201. AI(2,2)=AI(1,1)
  202. AI(2,3)=AI(1,1)
  203. AI(3,1)=AI(1,1)
  204. AI(3,2)=AI(1,1)
  205. AI(3,3)=AI(1,1)
  206. GOTO 75
  207. ENDIF
  208. C
  209. C---------- Cas du critere reduit a un point ----------------
  210. C
  211. IF (ABS(DPHI2).LE.10E-10) THEN
  212. IAPEX=2
  213. C WRITE(*,*)'IAPEX ds ACTI3 =',IAPEX
  214. C WRITE(*,*)'Dans l element',IBB
  215. C WRITE(*,*)'et au point d intégration',IGAU
  216. DO 52 I=1,NSTRSS
  217. DO 52 J=1,NSTRSS
  218. AI(I,J)=0.D0
  219. 52 CONTINUE
  220. AI(1,1)=1./3.
  221. AI(1,2)=AI(1,1)
  222. AI(1,2)=AI(1,1)
  223. AI(1,3)=AI(1,1)
  224. AI(2,1)=AI(1,1)
  225. AI(2,2)=AI(1,1)
  226. AI(2,3)=AI(1,1)
  227. AI(3,1)=AI(1,1)
  228. AI(3,2)=AI(1,1)
  229. AI(3,3)=AI(1,1)
  230. GOTO 75
  231. ENDIF
  232. C
  233. C ************** Mise a jour des contraintes ***************
  234. C
  235. C ---------------- calcul de la matrice A ------------------
  236. C
  237. DO 60 I=1,NSTRSS
  238. DO 60 J=1,NSTRSS
  239. A(I,J)=0.D0
  240. 60 CONTINUE
  241. C
  242. DG=YOUN/(1.D0+XNU)
  243. C
  244. A(1,1)=1.D0+2.*(DLAM1(1)*DG)/2.D0/DPHI1
  245. A(1,1)=A(1,1)+2.*(DLAM1(2)*DG)/2.D0/DPHI2
  246. A(2,2)=A(1,1)
  247. A(3,3)=A(1,1)
  248. A(1,2)=-(DLAM1(1)*DG)/2.D0/DPHI1-(DLAM1(2)*DG)/2.D0/DPHI2
  249. A(1,3)=A(1,2)
  250. A(2,1)=A(1,2)
  251. A(2,3)=A(1,2)
  252. A(3,1)=A(1,2)
  253. A(3,2)=A(1,2)
  254. A(4,4)=1.D0+3.*(DLAM1(1)*DG)/2.D0/DPHI1
  255. A(4,4)=A(4,4)+3.*(DLAM1(2)*DG)/2.D0/DPHI2
  256. A(5,5)=A(4,4)
  257. A(6,6)=A(4,4)
  258. C
  259. C -------------- invertion de la matrice A -----------------
  260. C
  261. CALL ZERO(AI,6,6)
  262. CALL ZERO(AIM,6,6)
  263. C
  264. DO 70 I=1,3
  265. DO 70 J=1,3
  266. AIM(I,J)=A(I,J)
  267. 70 CONTINUE
  268. CALL INVMA2(AIM,3,ISING)
  269. IF (ISING.EQ.1) THEN
  270. WRITE(*,*)'MATRICE AIM singuliere ds ACTI3'
  271. ENDIF
  272. DO 72 I=1,3
  273. DO 72 J=1,3
  274. AI(I,J)=AIM(I,J)
  275. 72 CONTINUE
  276. AI(4,4) = 1./A(4,4)
  277. AI(5,5) = 1./A(5,5)
  278. AI(6,6) = 1./A(6,6)
  279. C
  280. C -------------- mise a jour des contraintes ------------
  281. C
  282. 75 CONTINUE
  283. C
  284. DO 80 I=1,NSTRSS
  285. DEPSI(I)=SIGE(I)-DLAM1(1)*VEC1(I)
  286. DEPSI(I)=DEPSI(I)-DLAM1(2)*VEC2(I)
  287. 80 CONTINUE
  288. C
  289. DO 90 I=1,NSTRSS
  290. SIGF(I)=0.0D+00
  291. DO 90 J=1,NSTRSS
  292. SIGF(I)=SIGF(I)+AI(I,J)*DEPSI(J)
  293. 90 CONTINUE
  294. C
  295. C ******** Verification des criteres ****************
  296. C
  297. CALL DRUTRA(SIGF,SEQTT,BETINSA)
  298. FCRI1(1) = SEQTT - SEQ1
  299. CALL DRUCOM(SIGF,SEQCC,BETINSA)
  300. FCRI1(2) = SEQCC - SEQ2
  301. C
  302. IF (IAPEX.EQ.2) FCRI1(1)=0.D0
  303. C
  304. IF (IBROY.EQ.0.AND.(ABS(FCRI1(1)).GE.CRIMAX.OR
  305. * .ABS(FCRI1(2)).GE.CRIMAX)) THEN
  306. C WRITE(*,*)'****************************************'
  307. C WRITE(*,*)'LE RESIDU DIVERGE AVEC BROYDEN'
  308. C WRITE(*,*)'on passe donc a la secante'
  309. C WRITE(*,*)'Dans l element',IBB
  310. C WRITE(*,*)'et au point d intégration',IGAU
  311. C WRITE(*,*)'CRIMAX=',CRIMAX
  312. C WRITE(*,*)'****************************************'
  313. ITER=ITR
  314. ENDIF
  315. C
  316. C ******* Compteur sur la methode de resolution ****
  317. C
  318. IF (IBROY.EQ.0.AND.ITER.EQ.ITR) THEN
  319. IBROY=1
  320. ITANG=1
  321. ITER=1
  322. IAPEX=0
  323. IET=0
  324. IEC=0
  325. GOTO 12
  326. ENDIF
  327. C
  328. C ******* non convergence **************************
  329. C
  330. IF ((ABS(FCRI1(1)).GT.PRB.OR.ABS(FCRI1(2)).GT.PRB)
  331. *.AND.ITER.LT.ITR) THEN
  332. IF (IBROY.EQ.0) THEN
  333. CALL BROYDI(DJI,DLAM0,DLAM1,FCRI0,FCRI1)
  334. DLAM0(1)=DLAM1(1)
  335. DLAM0(2)=DLAM1(2)
  336. FCRI0(1)=FCRI1(1)
  337. FCRI0(2)=FCRI1(2)
  338. IF (ITER.GE.(ITR-1)) THEN
  339. C WRITE(*,*)'***********************'
  340. C WRITE(*,*)'BROYDEN n a pas aboutit'
  341. C WRITE(*,*)'ITER=',ITER
  342. C WRITE(*,*)'FCRIT=',FCRI0(1)
  343. C WRITE(*,*)'FCRIC=',FCRI0(2)
  344. C WRITE(*,*)'Dans l element',IBB
  345. C WRITE(*,*)'et au point d intégration',IGAU
  346. C WRITE(*,*)'***********************'
  347. ENDIF
  348. ITER=ITER+1
  349. GOTO 40
  350. ENDIF
  351. IF (IBROY.EQ.1.AND.ITANG.EQ.1) THEN
  352. DLAM0(1)=DLAM1(1)
  353. DLAM0(2)=DLAM1(2)
  354. FCRI0(1)=FCRI1(1)
  355. FCRI0(2)=FCRI1(2)
  356. IF (ITER.GE.(ITR-5)) THEN
  357. C WRITE(*,*)'ITER=',ITER
  358. C WRITE(*,*)'FCRIT=',FCRI0(1)
  359. C WRITE(*,*)'FCRIC=',FCRI0(2)
  360. ENDIF
  361. ITER=ITER+1
  362. GOTO 18
  363. ENDIF
  364. ENDIF
  365. IF (ITER.GE.ITR.AND.(ABS(FCRI1(1)).GT.PRB2.OR.
  366. * ABS(FCRI1(2)).GT.PRB2)) THEN
  367. WRITE(*,*)'NON CONVERGENCE INTERNE dans BEHAV3'
  368. WRITE(*,*)'Dans l element',IBB
  369. WRITE(*,*)'et au point d intégration',IGAU
  370. WRITE(*,*)'FCRIT=',FCRI0(1)
  371. WRITE(*,*)'FCRIC=',FCRI0(2)
  372. WRITE(*,*)'IPLA=',IPLA
  373. WRITE(*,*)'IFIS=',IFIS
  374. C STOP
  375. ENDIF
  376. C
  377. C **************** Fin des iterations internes *******************
  378. C
  379. DKT=DK1
  380. DKC=DK2
  381. SEQT=SEQ1
  382. SEQC=SEQ2
  383. C
  384. C ********************************************************************
  385. 100 CONTINUE
  386. IET=0
  387. IEC=0
  388. CRIMAX=0.D0
  389. C
  390. RETURN
  391. END
  392.  
  393.  
  394.  
  395.  
  396.  
  397.  
  398.  
  399.  

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