Télécharger betokt.eso

Retour à la liste

Numérotation des lignes :

  1. C BETOKT SOURCE PV 11/03/07 21:15:10 6885
  2. SUBROUTINE BETOKT(XMAT,TRAC,VARF,SIGF,NCURVT,NCURVC,XKTAT,KERRE)
  3. C
  4. C Matrice tangente pour le modele endommageable en cyclique
  5. C D. Combescure 18/09/95 ELSA-Ispra
  6. C
  7. C ENTREES:
  8. C -------
  9. C XMAT = Caracteristiques materiaux
  10. C
  11. C VARF = Variables internes
  12. C
  13. C SIGF = Contraintes
  14. C
  15. C SORTIES:
  16. C -------
  17. C XKTAT = Matrice tangente
  18. C
  19. C================================================================
  20. IMPLICIT INTEGER(I-N)
  21. IMPLICIT REAL*8(A-H,O-Z)
  22. -INC CCOPTIO
  23. DIMENSION XMAT(*)
  24. DIMENSION TRAC(*)
  25. C
  26. DIMENSION VARF(14),SIGF(4),EPSI(2),FG(2)
  27. DIMENSION XNOR1(4),XNOR2(4),XCOU1(3,3),XCOU2(3,3)
  28. DIMENSION GAMBDA(2),EEE(2,2),SIGP(2),EEE2(2,2)
  29. DIMENSION UUU(3,2),VVV(2,3),UUU2(3,2),VVV2(2,3),UUU3(3,2)
  30. DIMENSION HHH1(3,3),HHH2(3,3),HHH3(3,3)
  31. DIMENSION XKTAT(4,4),RESIST(6),XKISO(2)
  32. PARAMETER (UN=1.D0,UNDEMI=0.5D0,EPSILO=1.D-20,XZER=0.D0)
  33. PARAMETER (PRECIS=1.D-4)
  34. KERRE=0
  35. C================================================================
  36. C LECTURE DES CARACTERISTIQUES MAT. ET DES VAR. INT.
  37. C===============================================================
  38. C Pour la lecture des courbes
  39. IPT = 1
  40. IPC = IPT+2*NCURVT
  41. C
  42. YOUN = XMAT(1)
  43. XNU = XMAT(2)
  44. XGGG = (UNDEMI*YOUN/(UN+XNU))
  45. XHH1 = XMAT(3)
  46. DEGRAD = XMAT(8)/XMAT(9)
  47. RESIST(1) = XMAT(4)
  48. RESIST(2) = XMAT(4)
  49. RESIST(3) = XMAT(5)
  50. RESIST(4) = XMAT(5)
  51. RESIST(5) = XMAT(6)
  52. RESIST(6) = XMAT(7)
  53. XLCAT = XMAT(10)
  54. XLCAC = XMAT(11)
  55. CPHI = VARF(5) + UN
  56. SPHI = VARF(6)
  57. EPSI(1) = VARF(3)
  58. EPSI(2) = VARF(4)
  59. ICAS1 = nint(VARF(11))
  60. ICAS2 = nint(VARF(12))
  61. C==============================================================
  62. C RAIDEUR ELASTIQUE
  63. C===============================================================
  64. CALL ZERO(XKTAT,4,4)
  65. XKTAT(1,1) = YOUN/(UN - (XNU**2))
  66. XKTAT(1,2) = XNU*XKTAT(1,1)
  67. XKTAT(2,1) = XKTAT(1,2)
  68. XKTAT(2,2) = XKTAT(1,1)
  69. XKTAT(4,4) = YOUN/(2.D0*(UN+XNU))
  70. C==============================================================
  71. IF (ICAS1.EQ.0) THEN
  72. GOTO 9999
  73. ENDIF
  74. C
  75. CALL VALEUR(EPSI(1),TRAC(IPT),NCURVT,XLCAT,
  76. &XKISO(1),XPEN1,KERRE)
  77. CALL VALEUR(EPSI(2),TRAC(IPC),NCURVC,XLCAC,
  78. &XKISO(2),XPEN2,KERRE)
  79. GAMBDA(1) = VARF(13)
  80. GAMBDA(2) = VARF(14)
  81. IF (ICAS1.LE.4) THEN
  82. FG(1) = RESIST(ICAS1)
  83. ELSE
  84. FG(1) = XKISO(ICAS1-4)*RESIST(ICAS1)
  85. ENDIF
  86. IF (ICAS2.LE.4) THEN
  87. FG(2) = RESIST(ICAS2)
  88. ELSE
  89. FG(2) = XKISO(ICAS2-4)*RESIST(ICAS2)
  90. ENDIF
  91. C==============================================================
  92. C CALCUL DE L'INVERSE
  93. C==============================================================
  94. CALL ZERO(HHH1,3,3)
  95. HHH1(1,1) = UN/YOUN
  96. HHH1(1,2) = -XNU*HHH1(1,1)
  97. HHH1(2,1) = HHH1(1,2)
  98. HHH1(2,2) = HHH1(1,1)
  99. HHH1(3,3) = (2.D0*(UN+XNU)/YOUN)
  100. C
  101. CALL CONPRI (SIGF,CPHI,SPHI,SIGP,CPHI,SPHI)
  102. C
  103. IF (ICAS1.EQ.1) THEN
  104. ISGNTC = 1
  105. ELSE IF (ICAS1.EQ.2) THEN
  106. ISGNTC = 1
  107. ELSE IF (ICAS1.EQ.3) THEN
  108. ISGNTC = -1
  109. ELSE IF (ICAS1.EQ.4) THEN
  110. ISGNTC = -1
  111. ELSE IF (ICAS1.EQ.5) THEN
  112. ISGNTC = 1
  113. ELSE IF (ICAS1.EQ.6) THEN
  114. ISGNTC = -1
  115. ENDIF
  116. IF (ICAS1.LE.4) THEN
  117. IF (SIGP(2).GE.SIGP(1)) THEN
  118. IF ((ICAS1.EQ.1).OR.(ICAS1.EQ.3)) THEN
  119. ISGN2 = -1
  120. ELSE
  121. ISGN2 = 1
  122. ENDIF
  123. ELSE
  124. IF ((ICAS1.EQ.1).OR.(ICAS1.EQ.3)) THEN
  125. ISGN2 = 1
  126. ELSE
  127. ISGN2 = -1
  128. ENDIF
  129. ENDIF
  130. ELSE
  131. ISGN2 = 1
  132. ENDIF
  133. IF (ICAS1.LE.4) THEN
  134. CALL COURBU(SIGF,FG(1),ISGNTC,ISGN2,1,XCOU1)
  135. CALL NORMBL(SIGF,FG(1),ISGNTC,ISGN2,1,XNOR1)
  136. ELSE
  137. PHI = SQRT ((UNDEMI*(SIGF(1)-SIGF(2)))**2+SIGF(4)**2)
  138. IF ((PHI.LT.ABS(PRECIS*FG(1))).AND.(ICAS1.EQ.ICAS2)) THEN
  139. C Cas de l'apex ou on calcule directement H-1
  140. HHH1(1,1) = XGGG*((4.D0*PHI)+(YOUN*GAMBDA(1)))/
  141. & (2.D0*(UN-XNU)*((2.D0*PHI)+(XGGG*GAMBDA(1))))
  142. HHH1(2,2)=HHH1(1,1)
  143. HHH1(1,2) = XGGG*((4.D0*XNU*PHI)+(YOUN*GAMBDA(1)))/
  144. & (2.D0*(UN-XNU)*((2.D0*PHI)+(XGGG*GAMBDA(1))))
  145. HHH1(2,1)=HHH1(1,2)
  146. HHH1(3,3) = XGGG*PHI/(PHI+XGGG*GAMBDA(1))
  147. CALL NORMBL(SIGF,FG(1),ISGNTC,1,2,XNOR1)
  148. GOTO 1010
  149. ELSE
  150. CALL COURBU(SIGF,FG(1),ISGNTC,1,2,XCOU1)
  151. CALL NORMBL(SIGF,FG(1),ISGNTC,1,2,XNOR1)
  152. ENDIF
  153. ENDIF
  154. C
  155. C
  156. IF (ICAS1.EQ.ICAS2) THEN
  157. GOTO 1000
  158. ELSE
  159. IF (ICAS2.LE.4) THEN
  160. FG(2) = RESIST(ICAS2)
  161. ELSE
  162. FG(2) = XKISO(ICAS2-4)*RESIST(ICAS2)
  163. ENDIF
  164. IF (ICAS2.EQ.1) THEN
  165. ISGNTCB = 1
  166. ELSE IF (ICAS2.EQ.2) THEN
  167. ISGNTCB = 1
  168. ELSE IF (ICAS2.EQ.3) THEN
  169. ISGNTCB = -1
  170. ELSE IF (ICAS2.EQ.4) THEN
  171. ISGNTCB = -1
  172. ELSE IF (ICAS2.EQ.5) THEN
  173. ISGNTCB = 1
  174. ELSE IF (ICAS2.EQ.6) THEN
  175. ISGNTCB = -1
  176. ENDIF
  177. IF (ICAS2.LE.4) THEN
  178. IF (SIGP(2).GE.SIGP(1)) THEN
  179. IF ((ICAS2.EQ.1).OR.(ICAS2.EQ.3)) THEN
  180. ISGN2B = -1
  181. ELSE
  182. ISGN2B = 1
  183. ENDIF
  184. ELSE
  185. IF ((ICAS2.EQ.1).OR.(ICAS2.EQ.3)) THEN
  186. ISGN2B = 1
  187. ELSE
  188. ISGN2B = -1
  189. ENDIF
  190. ENDIF
  191. ELSE
  192. ISGN2B = 1
  193. ENDIF
  194. C
  195. IF (ICAS2.LE.4) THEN
  196. CALL COURBU(SIGF,FG(2),ISGNTCB,ISGN2B,1,XCOU2)
  197. CALL NORMBL(SIGF,FG(2),ISGNTCB,ISGN2B,1,XNOR2)
  198. ELSE
  199. CALL COURBU(SIGF,FG(2),ISGNTCB,1,2,XCOU2)
  200. CALL NORMBL(SIGF,FG(2),ISGNTCB,1,2,XNOR2)
  201. ENDIF
  202. C
  203. GOTO 2000
  204. ENDIF
  205. C
  206. C
  207. C
  208. 1000 CONTINUE
  209. C
  210. C Cas 1 mecanisme
  211. C
  212. C Calcul de l'inverse de H et de H
  213. CALL PLMATR(HHH1,UN,XCOU1,GAMBDA(1),3,3,HHH1)
  214. CALL INVALM(HHH1,3,3,IRD,1.D-8)
  215. IF (IRD.NE.0) THEN
  216. WRITE (IOIMP,*)'BETOKT: Determinant nul'
  217. GOTO 9999
  218. ENDIF
  219. C Point d'entree dans le cas de l'apex
  220. 1010 CONTINUE
  221. C Calcul de U,V et E-1
  222. CALL ZERO(UUU,3,2)
  223. CALL ZERO(VVV,2,3)
  224. CALL ZERO(UUU2,3,2)
  225. CALL ZERO(VVV2,2,3)
  226. UUU(1,1) = XNOR1(1)
  227. UUU(2,1) = XNOR1(2)
  228. UUU(3,1) = XNOR1(4)
  229. VVV(1,1) = XNOR1(1)
  230. VVV(1,2) = XNOR1(2)
  231. VVV(1,3) = XNOR1(4)
  232. C
  233. IF (ICAS1.LE.4) THEN
  234. EEE(1,1) = XHH1
  235. ELSE IF (ICAS1.EQ.5) THEN
  236. EEE(1,1) = ISGNTC*XPEN1*RESIST(ICAS1)
  237. ELSE IF (ICAS1.EQ.6) THEN
  238. EEE(1,1) = ISGNTC*XPEN2*RESIST(ICAS1)
  239. ENDIF
  240. C
  241. CALL MATMAT(HHH1,UUU,3,3,2,UUU2)
  242. CALL MATMAT(VVV,HHH1,2,3,3,VVV2)
  243. C
  244. CALL MATMAT(VVV,UUU2,2,3,2,HHH2)
  245. C
  246. EEE(1,1) = EEE(1,1) + HHH2(1,1)
  247. IF (ABS(EEE(1,1)).GT.EPSILO) THEN
  248. EEE(1,1) = -UN/EEE(1,1)
  249. ELSE
  250. WRITE (IOIMP,*)'BETOKT: Determinant nul'
  251. GOTO 9999
  252. ENDIF
  253. CALL MATMAT(UUU2,VVV2,3,2,3,HHH3)
  254. CALL PLMATR(HHH1,UN,HHH3,EEE(1,1),3,3,HHH1)
  255. C
  256. GOTO 3000
  257. C
  258. 2000 CONTINUE
  259. C
  260. C Cas 2 mecanismes
  261. C
  262. C Calcul de H-1 et H
  263. CALL PLMATR(HHH1,UN,XCOU1,GAMBDA(1),3,3,HHH1)
  264. CALL PLMATR(HHH1,UN,XCOU2,GAMBDA(2),3,3,HHH1)
  265. CALL INVALM(HHH1,3,3,IRD,1.D-8)
  266. IF (IRD.NE.0) THEN
  267. WRITE (IOIMP,*)'BETOKT: Determinant nul'
  268. GOTO 9999
  269. ENDIF
  270. C Calcul de U ,V et E
  271. UUU(1,1) = XNOR1(1)
  272. UUU(2,1) = XNOR1(2)
  273. UUU(3,1) = XNOR1(4)
  274. UUU(1,2) = XNOR2(1)
  275. UUU(2,2) = XNOR2(2)
  276. UUU(3,2) = XNOR2(4)
  277. VVV(1,1) = XNOR1(1)
  278. VVV(1,2) = XNOR1(2)
  279. VVV(1,3) = XNOR1(4)
  280. VVV(2,1) = XNOR2(1)
  281. VVV(2,2) = XNOR2(2)
  282. VVV(2,3) = XNOR2(4)
  283. C
  284. IF (ICAS1.LE.4) THEN
  285. EEE(1,1) = XHH1
  286. IF (ICAS2.EQ.6) THEN
  287. EEE(2,1) = DEGRAD*ISGNTCB*RESIST(6)
  288. &*ISGNTC*RESIST(ICAS1)*XPEN2
  289. ELSE
  290. EEE(2,1) = XZER
  291. ENDIF
  292. ELSE IF (ICAS1.EQ.5) THEN
  293. EEE(1,1) = ISGNTC*XPEN1*RESIST(ICAS1)
  294. EEE(2,1) = XZER
  295. ELSE IF (ICAS1.EQ.6) THEN
  296. EEE(1,1) = ISGNTC*XPEN2*RESIST(ICAS1)
  297. EEE(2,1) = XZER
  298. ENDIF
  299. IF (ICAS2.LE.4) THEN
  300. EEE(2,2) = XHH1
  301. IF (ICAS1.EQ.6) THEN
  302. EEE(1,2) = DEGRAD*ISGNTC*RESIST(6)
  303. &*ISGNTCB*RESIST(ICAS2)*XPEN2
  304. ELSE
  305. EEE(1,2) = XZER
  306. ENDIF
  307. ELSE IF (ICAS2.EQ.5) THEN
  308. EEE(2,2) = ISGNTCB*XPEN1*RESIST(ICAS2)
  309. EEE(1,2) = XZER
  310. ELSE IF (ICAS2.EQ.6) THEN
  311. EEE(2,2) = ISGNTCB*XPEN2*RESIST(ICAS2)
  312. EEE(1,2) = XZER
  313. ENDIF
  314. C
  315. CALL MATMAT(HHH1,UUU,3,3,2,UUU2)
  316. CALL MATMAT(VVV,HHH1,2,3,3,VVV2)
  317. C
  318. CALL MATMAT(VVV,UUU2,2,3,2,HHH2)
  319. C
  320. CALL PLMATR(EEE,UN,HHH2,UN,2,2,EEE)
  321. XDET = EEE(1,1)*EEE(2,2) - EEE(1,2)*EEE(2,1)
  322. IF (ABS(XDET).GT.EPSILO) THEN
  323. XDET = UN/XDET
  324. EEE2(1,1) = XDET*EEE(2,2)
  325. EEE2(2,2) = XDET*EEE(1,1)
  326. EEE2(2,1) = -XDET*EEE(1,2)
  327. EEE2(1,2) = -XDET*EEE(2,1)
  328. ELSE
  329. WRITE (IOIMP,*)'BETOKT: Determinant nul'
  330. GOTO 9999
  331. ENDIF
  332. CALL MATMAT(UUU2,EEE2,3,2,2,UUU3)
  333. CALL MATMAT(UUU3,VVV2,3,2,3,HHH3)
  334. CALL PLMATR(HHH1,UN,HHH3,-UN,3,3,HHH1)
  335. C
  336. GOTO 3000
  337. C
  338. C==============================================================
  339. 3000 CONTINUE
  340. CALL ZERO(XKTAT,4,4)
  341. XKTAT(1,1) = HHH1(1,1)
  342. XKTAT(1,2) = HHH1(1,2)
  343. XKTAT(2,1) = HHH1(2,1)
  344. XKTAT(2,2) = HHH1(2,2)
  345. XKTAT(4,4) = HHH1(3,3)
  346. XKTAT(4,1) = HHH1(3,1)
  347. XKTAT(4,2) = HHH1(3,2)
  348. XKTAT(1,4) = HHH1(1,3)
  349. XKTAT(2,4) = HHH1(2,3)
  350. 9999 CONTINUE
  351. C===============================================================
  352. C Fin de la routine
  353. C===============================================================
  354. RETURN
  355. END
  356.  
  357.  
  358.  
  359.  

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