Télécharger dohot1.eso

Retour à la liste

Numérotation des lignes :

  1. C DOHOT1 SOURCE BP208322 17/03/01 21:17:03 9325
  2.  
  3. SUBROUTINE DOHOT1(IVAMAT,NMATT,IVACON,NSTRS,IVARI,NVART,TRAC,
  4. & LTRAC,IGAU,IB,MATE,MAPL,XPREC,DTPS,IFOU,LHOOK,
  5. & DDHOOK,IRET)
  6. *_______________________________________________________________________
  7. *
  8. * entrees :
  9. * --------
  10. *
  11. * ivamat pointeur sur un segment mptval contenant les materiaux
  12. * nmatt nombre de composantes du materiau
  13. * ivacon pointeur sur un segment mptval contenant les contraintes
  14. * nstrs nombre de composantes de contraintes
  15. * ivari pointeur sur un segment mptval contenant les variables
  16. * internes
  17. * nvart nombre de composantes de variables internes
  18. * trac courbe de traction
  19. * ltrac taille du tableau trac
  20. * igau numero du point de gauss ou l'on veut le mat. de hooke
  21. * ib numero de l'element ou l'on veut le mat. de hooke
  22. * mate numeror du materiau lineaire
  23. * mapl numero du materiau plastique
  24. * xprec flottant precision
  25. * dtps flottant pas de temps pour les modeles visqueux
  26. * ifou option ifour de ccoptio
  27. * lhook taille de la matrice de hooke
  28. *
  29. * sorties :
  30. * --------
  31. *
  32. * ddhook matrice de hooke tangente pour l'element ib
  33. * iret = 1 si ok
  34. * = < 0 : no de l erreur a appeler dans ktanga
  35. * = 0 : autre erreur
  36. *
  37. * passage aux nouveaux chamelem par jm campenon le 05/91
  38. *_______________________________________________________________________
  39. *
  40. IMPLICIT INTEGER(I-N)
  41. IMPLICIT REAL*8(A-H,O-Z)
  42. *
  43. -INC CCOPTIO
  44. -INC SMCHAML
  45. *
  46. SEGMENT MPTVAL
  47. INTEGER IPOS(NS) ,NSOF(NS)
  48. INTEGER IVAL(NCOSOU)
  49. CHARACTER*16 TYVAL(NCOSOU)
  50. ENDSEGMENT
  51. *
  52. SEGMENT WRK5
  53. INTEGER NTRAT,NTRAC
  54. ENDSEGMENT
  55. *
  56. SEGMENT WTRAV
  57. REAL*8 TXR(IDIM,IDIM)
  58. ENDSEGMENT
  59.  
  60. SEGMENT WRKMAT
  61. REAL*8 VALMAT(NMATT)
  62. REAL*8 D1HO(LHOOK,LHOOK),ROTH(LHOOK,LHOOK)
  63. REAL*8 XLOC(3,3),XGLOB(3,3)
  64. ENDSEGMENT
  65. *
  66. PARAMETER(XZER=0.D0,UN=1.D0,DEUX=2.D0,TROIS=3.D0,QUATRE=4.D0)
  67. PARAMETER(CINQ=5.D0,SIX=6.D0,HUIT=8.D0,XNEUF=9.D0,DOUZE=12.D0)
  68. PARAMETER(UNTIER=.33333 33333 33333D0,TRDEMI=1.5D0)
  69. PARAMETER(DETIER=.66666 66666 66666D0,UNDEMI=0.5D0)
  70. PARAMETER(UNQU=.25D0)
  71. *
  72. DIMENSION DDHOOK(LHOOK,*)
  73. DIMENSION TRAC(*)
  74. *
  75. * arguments pour l appel a DOHMAS et DOHMAO
  76. DIMENSION VELA(2)
  77. CHARACTER*8 CMATE
  78. *
  79. DIMENSION XSTRS(8),DEVIAT(8),XMAT(30)
  80. DIMENSION DDF(8),DDG(8),DHDF(8),DHDG(8)
  81. DIMENSION XVAR(14),XKTA(3,3)
  82. DIMENSION PMATRIC(6,6),GMATRIC(6,6),GDMATR(6,6),GINV(6,6),COEF3(6)
  83. * matrice P pour le modele viscoplastique parfait
  84. * la matrice sert a passer de la contrainte a la contrainte deviatorique
  85. DATA PMATRIC /2.D0,-1.D0,-1.D0,0.D0,0.D0,0.D0,
  86. & -1.D0,2.D0,-1.D0,0.D0,0.D0,0.D0,
  87. & -1.D0,-1.D0,2.D0,0.D0,0.D0,0.D0,
  88. & 0.D0, 0.D0,0.D0,6.D0,0.D0,0.D0,
  89. & 0.D0, 0.D0,0.D0,0.D0,6.D0,0.D0,
  90. & 0.D0, 0.D0,0.D0,0.D0,0.D0,6.D0/
  91. *
  92. DATA COEF3 /1.D0,1.D0,1.D0,2.D0,2.D0,2.D0/
  93. *
  94. *-----------------------------------------------------------------------
  95. * formulation massive (MFR=1)
  96. * Pour autres formulations (ex : coques) voir dans dohota.eso
  97. *-----------------------------------------------------------------------
  98. IF ((MATE.EQ.1).AND.(MAPL.EQ.54)) THEN
  99. WRK5 = IRET
  100. ELSE IF ((MATE.EQ.2).OR.(MATE.EQ.3).OR.(MATE.EQ.4)) THEN
  101. WTRAV = IRET
  102. ENDIF
  103. *
  104. IRET=0
  105. CALL ZERO(DDHOOK,LHOOK,LHOOK)
  106. CALL ZERO(XMAT,30,1)
  107. *
  108. *-----------------------------------------------------------------------
  109. * MATERIAU ISOTROPE
  110. *-----------------------------------------------------------------------
  111. IF (MATE.EQ.1) THEN
  112. *-----------------------------------------------------------------------
  113. C Initialisations
  114. XXHH =0.D0
  115. SIGY =0.D0
  116. ALP =0.D0
  117. C Initialisation de ILOGEL : sortie de PENT1
  118. C =1 on doit/peut calculer DDHOOK tangent ; =0 DDHOOK elastique
  119. ILOGEL = 0
  120. *
  121. C Contraintes
  122. MPTVAL=IVACON
  123. DO 50 ICOMP=1,NSTRS
  124. MELVAL=IVAL(ICOMP)
  125. IGMN=MIN(IGAU,VELCHE(/1))
  126. IBMN=MIN(IB ,VELCHE(/2))
  127. XSTRS(ICOMP)=VELCHE(IGMN,IBMN)
  128. 50 CONTINUE
  129. C
  130. C Deformation plastique
  131. C
  132. MPTVAL=IVARI
  133. MELVAL=IVAL(1)
  134. IGMN=MIN(IGAU,VELCHE(/1))
  135. IBMN=MIN(IB ,VELCHE(/2))
  136. EPS =VELCHE(IGMN,IBMN)
  137. C si autres variables internes
  138. IF (NVART.EQ.NSTRS+1) THEN
  139. DO 60 ICOMP=1,NSTRS
  140. MELVAL=IVAL(ICOMP+1)
  141. IGMN=MIN(IGAU,VELCHE(/1))
  142. IBMN=MIN(IB ,VELCHE(/2))
  143. XSTRS(ICOMP)= XSTRS(ICOMP)-VELCHE(IGMN,IBMN)
  144. 60 CONTINUE
  145. ENDIF
  146. C
  147. C Recuperation des caracteristiques elastiques isotropes
  148. C
  149. MPTVAL=IVAMAT
  150. C Module d Young
  151. MELVAL=IVAL(1)
  152. IGMN=MIN(IGAU,VELCHE(/1))
  153. IBMN=MIN(IB ,VELCHE(/2))
  154. YOU =VELCHE(IGMN,IBMN)
  155. C Coefficient de Poisson
  156. MELVAL=IVAL(2)
  157. IGMN=MIN(IGAU,VELCHE(/1))
  158. IBMN=MIN(IB ,VELCHE(/2))
  159. XNU =VELCHE(IGMN,IBMN)
  160. *
  161. *-----------------------------------------------------------------------
  162. * Recuperation des caracteristiques
  163. * Calcul eventuel de la contrainte equivalente
  164. *-----------------------------------------------------------------------
  165. C Plastique parfait
  166. IF (MAPL.EQ.1) THEN
  167. *
  168. MPTVAL=IVAMAT
  169. MELVAL=IVAL(3)
  170. IGMN=MIN(IGAU,VELCHE(/1))
  171. IBMN=MIN(IB ,VELCHE(/2))
  172. SIGY =VELCHE(IGMN,IBMN)
  173. *
  174. C Calcul de la contrainte equivalente
  175. XXXX = XSTRS(1)*XSTRS(1) + XSTRS(2)*XSTRS(2)
  176. & + XSTRS(3)*XSTRS(3)
  177. YYYY = XSTRS(1)*XSTRS(2) + XSTRS(2)*XSTRS(3)
  178. & + XSTRS(3)*XSTRS(1)
  179. SIGET2 = XXXX - YYYY
  180. IF (NSTRS.GE.4) THEN
  181. DO 61 I=4,NSTRS
  182. SIGET2 = SIGET2 + TROIS*XSTRS(I)*XSTRS(I)
  183. 61 CONTINUE
  184. ENDIF
  185. SIGET = SQRT(SIGET2)
  186. *
  187. C A t on plastifie ?
  188. CALL PENT1(EPS,SIGET,YOU,SIGY,XXHH,TRAC,LTRAC,MAPL,XPREC,
  189. & YOUTA,ILOGEL)
  190. IF (ILOGEL.EQ.-1) THEN
  191. IRET=-275
  192. RETURN
  193. ENDIF
  194. *
  195. C Drucker simple : ALFA * Tr(S) + Seq = K
  196. ELSE IF (MAPL.EQ.3) THEN
  197. *
  198. MPTVAL=IVAMAT
  199. MELVAL=IVAL(3)
  200. IGMN=MIN(IGAU,VELCHE(/1))
  201. IBMN=MIN(IB ,VELCHE(/2))
  202. XLTR =VELCHE(IGMN,IBMN)
  203. MELVAL=IVAL(4)
  204. IGMN=MIN(IGAU,VELCHE(/1))
  205. IBMN=MIN(IB ,VELCHE(/2))
  206. XLCS =VELCHE(IGMN,IBMN)
  207. *
  208. DEN = ABS(XLCS) + XLTR
  209. IF(DEN.EQ.0.D0) THEN
  210. IRET=-366
  211. RETURN
  212. ENDIF
  213. C ALP = ALFA ; XXHH = BETA = 1 ; SIGY = K
  214. ALP = (ABS(XLCS) - XLTR)/DEN
  215. XXHH = UN
  216. SIGY = DEUX*ABS(XLCS)*XLTR/DEN
  217. *
  218. * XMAT(1) = 2.0D0*ABS(ALP)*SIGY/DEN
  219. * XMAT(2) = (ABS(ALP) - SIGY)/DEN
  220. * XMAT(3) = 1.D0
  221. XMAT(1) = ALP
  222. XMAT(2) = XXHH
  223. XMAT(3) = SIGY
  224. XMAT(4)=XMAT(1)
  225. XMAT(5)=XMAT(2)
  226. XMAT(6)=XMAT(1)
  227. XMAT(7)=XMAT(2)
  228. XMAT(8)=XMAT(3)
  229. XMAT(9)=0.D0
  230. *
  231. C Calcul de la contrainte equivalente
  232. XXXX = XSTRS(1)*XSTRS(1) + XSTRS(2)*XSTRS(2)
  233. & + XSTRS(3)*XSTRS(3)
  234. YYYY = XSTRS(1)*XSTRS(2) + XSTRS(2)*XSTRS(3)
  235. & + XSTRS(3)*XSTRS(1)
  236. SIGET2 = XXXX - YYYY
  237. IF (NSTRS.GE.4) THEN
  238. DO 62 I=4,NSTRS
  239. SIGET2 = SIGET2 + TROIS*XSTRS(I)*XSTRS(I)
  240. 62 CONTINUE
  241. ENDIF
  242. SIGET = SQRT(SIGET2)
  243. *
  244. C A t on plastifie ?
  245. SIDD = ALP * (XSTRS(1)+XSTRS(2)+XSTRS(3)) + SIGET
  246. CALL PENT1(EPS,SIDD,YOU,SIGY,XXHH,TRAC,LTRAC,MAPL,XPREC,
  247. & YOUTA,ILOGEL)
  248. IF (ILOGEL.EQ.-1) THEN
  249. IRET=-275
  250. RETURN
  251. ENDIF
  252. *
  253. C Cinematique
  254. ELSE IF (MAPL.EQ.4) THEN
  255. *
  256. MPTVAL=IVAMAT
  257. MELVAL=IVAL(3)
  258. IGMN=MIN(IGAU,VELCHE(/1))
  259. IBMN=MIN(IB ,VELCHE(/2))
  260. SIGY =VELCHE(IGMN,IBMN)
  261. *
  262. MELVAL=IVAL(4)
  263. IGMN=MIN(IGAU,VELCHE(/1))
  264. IBMN=MIN(IB ,VELCHE(/2))
  265. XXHH=VELCHE(IGMN,IBMN)
  266. *
  267. C Calcul de la contrainte equivalente
  268. XXXX = XSTRS(1)*XSTRS(1) + XSTRS(2)*XSTRS(2)
  269. & + XSTRS(3)*XSTRS(3)
  270. YYYY = XSTRS(1)*XSTRS(2) + XSTRS(2)*XSTRS(3)
  271. & + XSTRS(3)*XSTRS(1)
  272. SIGET2 = XXXX - YYYY
  273. IF (NSTRS.GE.4) THEN
  274. DO 63 I=4,NSTRS
  275. SIGET2 = SIGET2 + TROIS*XSTRS(I)*XSTRS(I)
  276. 63 CONTINUE
  277. ENDIF
  278. SIGET = SQRT(SIGET2)
  279. *
  280. C A t on plastifie ?
  281. CALL PENT1(EPS,SIGET,YOU,SIGY,XXHH,TRAC,LTRAC,MAPL,XPREC,
  282. & YOUTA,ILOGEL)
  283. IF (ILOGEL.EQ.-1) THEN
  284. IRET=-275
  285. RETURN
  286. ENDIF
  287. *
  288. C Plastique isotrope
  289. ELSE IF (MAPL.EQ.5) THEN
  290. *
  291. C Calcul de la contrainte equivalente
  292. XXXX = XSTRS(1)*XSTRS(1) + XSTRS(2)*XSTRS(2)
  293. & + XSTRS(3)*XSTRS(3)
  294. YYYY = XSTRS(1)*XSTRS(2) + XSTRS(2)*XSTRS(3)
  295. & + XSTRS(3)*XSTRS(1)
  296. SIGET2 = XXXX - YYYY
  297. IF (NSTRS.GE.4) THEN
  298. DO 64 I=4,NSTRS
  299. SIGET2 = SIGET2 + TROIS*XSTRS(I)*XSTRS(I)
  300. 64 CONTINUE
  301. ENDIF
  302. SIGET = SQRT(SIGET2)
  303. *
  304. C A t on plastifie ?
  305. CALL PENT1(EPS,SIGET,YOU,SIGY,XXHH,TRAC,LTRAC,MAPL,XPREC,
  306. & YOUTA,ILOGEL)
  307. IF (ILOGEL.EQ.-1) THEN
  308. IRET=-275
  309. RETURN
  310. ENDIF
  311. *
  312. C Drucker Prager
  313. ELSE IF (MAPL.EQ.15) THEN
  314. *
  315. MPTVAL=IVAMAT
  316. DO 1106 IC=3,NMATT
  317. MELVAL=IVAL(IC)
  318. XMAT(IC-2)=VELCHE(IGMN,IBMN)
  319. 1106 CONTINUE
  320. *
  321. ** ALP = ALFA ou ETA ; XXHH = BETA ou MU ; SIGY = K ou KL+H*EPS
  322. IF(EPS.EQ.0.D0) THEN
  323. ALP=XMAT(6)
  324. XXHH=XMAT(7)
  325. SIGY=XMAT(8)
  326. ELSE
  327. ALP=XMAT(1)
  328. XXHH=XMAT(2)
  329. SIGY=XMAT(3)+XMAT(9)*EPS
  330. ENDIF
  331. *
  332. C Calcul de la contrainte equivalente
  333. XXXX = XSTRS(1)*XSTRS(1) + XSTRS(2)*XSTRS(2)
  334. & + XSTRS(3)*XSTRS(3)
  335. YYYY = XSTRS(1)*XSTRS(2) + XSTRS(2)*XSTRS(3)
  336. & + XSTRS(3)*XSTRS(1)
  337. SIGET2 = XXXX - YYYY
  338. IF (NSTRS.GE.4) THEN
  339. DO 66 I=4,NSTRS
  340. SIGET2 = SIGET2 + TROIS*XSTRS(I)*XSTRS(I)
  341. 66 CONTINUE
  342. ENDIF
  343. SIGET = SQRT(SIGET2)
  344. *
  345. C A t on plastifie ?
  346. SIDD = ALP * (XSTRS(1)+XSTRS(2)+XSTRS(3)) + XXHH * SIGET
  347. CALL PENT1(EPS,SIDD,YOU,SIGY,XXHH,TRAC,LTRAC,MAPL,XPREC,
  348. & YOUTA,ILOGEL)
  349. IF (ILOGEL.EQ.-1) THEN
  350. IRET=-275
  351. RETURN
  352. ENDIF
  353. *
  354. C Viscoplastique parfait
  355. ELSE IF (MAPL.EQ.43) THEN
  356. *
  357. C Recuperation des caracteristiques
  358. MPTVAL=IVAMAT
  359. MELVAL = IVAL(3)
  360. IGMN=MIN(IGAU,VELCHE(/1))
  361. IBMN=MIN(IB ,VELCHE(/2))
  362. SIGY = VELCHE(IGMN,IBMN)
  363. MELVAL = IVAL(5)
  364. IGMN=MIN(IGAU,VELCHE(/1))
  365. IBMN=MIN(IB ,VELCHE(/2))
  366. CK = VELCHE(IGMN,IBMN)
  367. MELVAL = IVAL(4)
  368. IGMN=MIN(IGAU,VELCHE(/1))
  369. IBMN=MIN(IB ,VELCHE(/2))
  370. CN = VELCHE(IGMN,IBMN)
  371. ILOGEL = 1
  372. *
  373. C Betocy
  374. ELSE IF (MAPL.EQ.54) THEN
  375. C
  376. MPTVAL=IVAMAT
  377. DO 1107 IC=1,NMATT-2
  378. MELVAL=IVAL(IC)
  379. IGMN=MIN(IGAU,VELCHE(/1))
  380. IBMN=MIN(IB ,VELCHE(/2))
  381. XMAT(IC)=VELCHE(IGMN,IBMN)
  382. 1107 CONTINUE
  383. C Lecture Variables internes
  384. MPTVAL=IVARI
  385. DO 1108 IC=1,NVART
  386. MELVAL=IVAL(IC)
  387. IGMN=MIN(IGAU,VELCHE(/1))
  388. IBMN=MIN(IB ,VELCHE(/2))
  389. XVAR(IC)=VELCHE(IGMN,IBMN)
  390. 1108 CONTINUE
  391. ILOGEL = 1
  392. *
  393. C Rotating Crack
  394. ELSE IF (MAPL.EQ.55) THEN
  395. C Lecture materiau
  396. MPTVAL=IVAMAT
  397. DO 1109 IC=1,NMATT
  398. MELVAL=IVAL(IC)
  399. IGMN=MIN(IGAU,VELCHE(/1))
  400. IBMN=MIN(IB ,VELCHE(/2))
  401. XMAT(IC)=VELCHE(IGMN,IBMN)
  402. 1109 CONTINUE
  403. C Lecture Variables internes
  404. MPTVAL=IVARI
  405. DO 1110 IC=1,NVART
  406. MELVAL=IVAL(IC)
  407. IGMN=MIN(IGAU,VELCHE(/1))
  408. IBMN=MIN(IB ,VELCHE(/2))
  409. XVAR(IC)=VELCHE(IGMN,IBMN)
  410. 1110 CONTINUE
  411. ILOGEL = 1
  412. *
  413. * BCN (start)
  414. ELSE IF ((MAPL.EQ.111).OR.(MAPL.EQ.112).OR.(MAPL.EQ.113)) THEN
  415. IF (MAPL.EQ.112.OR.MAPL.EQ.113) THEN
  416. C J2 (MAPL=112) or Rounded Hyperbolic Mohr-Coulomb (MAPL=113)
  417. C get auxiliary variable for CTM
  418. MPTVAL=IVARI
  419. MELVAL=IVAL(2)
  420. IGMN=MIN(IGAU,VELCHE(/1))
  421. IBMN=MIN(IB ,VELCHE(/2))
  422. AUX1 =VELCHE(IGMN,IBMN)
  423. ELSE IF (MAPL.EQ.111) THEN
  424. C MRS-Lade (MAPL=111)
  425. C get auxiliary variable for CTM
  426. MPTVAL=IVARI
  427. MELVAL=IVAL(1)
  428. IGMN=MIN(IGAU,VELCHE(/1))
  429. IBMN=MIN(IB ,VELCHE(/2))
  430. EPSCON =VELCHE(IGMN,IBMN)
  431. MPTVAL=IVARI
  432. MELVAL=IVAL(2)
  433. IGMN=MIN(IGAU,VELCHE(/1))
  434. IBMN=MIN(IB ,VELCHE(/2))
  435. EPSCAP =VELCHE(IGMN,IBMN)
  436. MPTVAL=IVARI
  437. MELVAL=IVAL(3)
  438. IGMN=MIN(IGAU,VELCHE(/1))
  439. IBMN=MIN(IB ,VELCHE(/2))
  440. AUX1CON =VELCHE(IGMN,IBMN)
  441. MPTVAL=IVARI
  442. MELVAL=IVAL(4)
  443. IGMN=MIN(IGAU,VELCHE(/1))
  444. IBMN=MIN(IB ,VELCHE(/2))
  445. AUX1CAP =VELCHE(IGMN,IBMN)
  446. ENDIF
  447. C The three models: read material constants in XMAT
  448. XMAT(1)=YOU
  449. XMAT(2)=XNU
  450. XMAT(3)=0.D0
  451. XMAT(4)=0.D0
  452. MPTVAL=IVAMAT
  453. DO IC=3,NMATT
  454. MELVAL=IVAL(IC)
  455. IGMN=MIN(IGAU,VELCHE(/1))
  456. IBMN=MIN(IB ,VELCHE(/2))
  457. XMAT(IC+2)=VELCHE(IGMN,IBMN)
  458. ENDDO
  459. ILOGEL = 1
  460. ENDIF
  461. * BCN (end)
  462. *
  463. *-----------------------------------------------------------------------
  464. * Calcul de la matrice elastique
  465. *-----------------------------------------------------------------------
  466. VELA(1) = YOU
  467. VELA(2) = XNU
  468. * MATE = 1 --> CMATE = ISOTROPE
  469. CMATE = 'ISOTROPE'
  470. CALL DOHMAS(VELA,CMATE,IFOU,LHOOK,1,DDHOOK,IRETOU)
  471. IF (IRETOU.EQ.0) THEN
  472. IRET = -81
  473. RETURN
  474. ENDIF
  475. IRET = 1
  476. *
  477. * ILOGEL = 0 on a fini, =1 on calcule la matrice tangente
  478. IF (ILOGEL.EQ.0) RETURN
  479. * series de Fourier, on rend la matrice elastique
  480. IF (IFOU.EQ.1) RETURN
  481. *
  482. *-----------------------------------------------------------------------
  483. * Calcul de la matrice tangente
  484. *-----------------------------------------------------------------------
  485. c MODELE MRS_LADE
  486. IF (MAPL.EQ.111) THEN
  487. *
  488. * notice de l operateur MODE: "option de calcul : plan defo,axis"
  489. * sinon, on rend la matrice elastique
  490. IF ((IFOU.EQ.-1).OR.(IFOU.EQ.0)) THEN
  491. nescri =0
  492. nues =6
  493. C mrs-lade requiere siempre derivacion numerica
  494. nnumer=3
  495. deltax=2.D0**(int(log10(1.D-6)/log10(2.D0)))
  496. call mtc_MRSMAC(DDHOOK,LHOOK,XSTRS,NSTRS,EPSCON,EPSCAP,
  497. & AUX1CON,AUX1CAP,XMAT,nescri,nues,
  498. & nnumer,deltax,kerre)
  499. * kerre = 1 : formulacion no disponible
  500. IF (KERRE.EQ.1) THEN
  501. IRET = -328
  502. ENDIF
  503. ENDIF
  504. *
  505. c MODELE J2
  506. ELSE IF (MAPL.EQ.112) THEN
  507. * notice de l operateur MODE: "option de calcul : plan defo,axis"
  508. * sinon, on rend la matrice elastique
  509. IF ((IFOU.EQ.-1).OR.(IFOU.EQ.0)) THEN
  510. nescri =0
  511. nues =6
  512. call mtc_j2(DDHOOK,LHOOK,XSTRS,NSTRS,EPS,AUX1,XMAT,
  513. & nescri,nues,kerre)
  514. * kerre vaut toujours 0
  515. ENDIF
  516. *
  517. c MODELE RH_COULOMB
  518. ELSE IF (MAPL.EQ.113) THEN
  519. * notice de l operateur MODE: "option de calcul : plan defo,axis"
  520. * sinon, on rend la matrice elastique
  521. IF ((IFOU.EQ.-1).OR.(IFOU.EQ.0)) THEN
  522. nescri =0
  523. nues =6
  524. call mtc_rhmc(DDHOOK,LHOOK,XSTRS,NSTRS,AUX1,XMAT,
  525. & nescri,nues,kerre)
  526. * kerre vaut toujours 0
  527. ENDIF
  528. *
  529. C Viscoplastique parfait
  530. ELSE IF (MAPL.EQ.43) THEN
  531. *
  532. * tout sauf contraintes planes
  533. * sinon on rend la matrice elastique
  534. IF ((IFOU.EQ.-1).OR.(IFOU.EQ.-3).OR.(IFOU.EQ.0).OR.
  535. & (IFOU.EQ.2).OR.(IFOU.EQ.3).OR.(IFOU.EQ.7).OR.(IFOU.EQ.9)
  536. & .OR.(IFOU.EQ.11).OR.(IFOU.EQ.12).OR.(IFOU.EQ.14).OR.
  537. & (IFOU.EQ.15)) THEN
  538. * SIGMA EQUIVALENT: (3/2(s':s'))**0.5D0
  539. XXXX = XSTRS(1)*XSTRS(1) + XSTRS(2)*XSTRS(2)
  540. & + XSTRS(3)*XSTRS(3)
  541. YYYY = XSTRS(1)*XSTRS(2) + XSTRS(2)*XSTRS(3)
  542. & + XSTRS(3)*XSTRS(1)
  543. SIGET2 = XXXX - YYYY
  544. IF (NSTRS.GE.4) THEN
  545. DO 67 I=4,NSTRS
  546. SIGET2 = SIGET2 + TROIS*XSTRS(I)*XSTRS(I)
  547. 67 CONTINUE
  548. ENDIF
  549. SIGEQ = SQRT(SIGET2)
  550. *
  551. XMOY=(XSTRS(1)+XSTRS(2)+XSTRS(3))*UNTIER
  552. DO 200 IA=1,3
  553. DEVIAT(IA)=XSTRS(IA)- XMOY
  554. 200 CONTINUE
  555. IF (NSTRS.GE.4) THEN
  556. DO 300 IA=4,NSTRS
  557. DEVIAT(IA)=XSTRS(IA)
  558. 300 CONTINUE
  559. ENDIF
  560. *
  561. SYN = ( SIGEQ - SIGY ) / CK
  562. IF (SYN .GT. XPREC) THEN
  563. * on ne travaille que si l'evolution est plastique
  564. *
  565. AUX6 = UNDEMI * YOU / (UN+XNU)
  566. *
  567. * ON CALCULE LA MATRICE G
  568. COEF1 = SYN ** CN / SIGEQ*DTPS*2.D0*AUX6/3.D0
  569. COEF2=SYN**(CN-1.D0)*(SYN/SIGEQ+CN/CK)/SIGET2*DTPS
  570. & *2.D0*AUX6
  571. DO 610 I=1,NSTRS
  572. DO 600 J=1,NSTRS
  573. GMATRIC(I,J)=COEF1*PMATRIC(I,J)+
  574. & COEF2*DEVIAT(I)*DEVIAT(J)*COEF3(I)
  575. 600 CONTINUE
  576. GMATRIC(I,I)=GMATRIC(I,I)+1.D0
  577. 610 CONTINUE
  578.  
  579. CALL INVHOK(GMATRIC,GINV,NSTRS)
  580. *
  581. * on calcule InvG . D
  582. DO 650 I=1,NSTRS
  583. DO 640 J=I,NSTRS
  584. GDMATR(I,J)=0.D0
  585. DO 630 K=1,NSTRS
  586. GDMATR(I,J)=GDMATR(I,J)+GINV(I,K)*DDHOOK(K,J)
  587. 630 CONTINUE
  588. GDMATR(J,I)=GDMATR(I,J)
  589. 640 CONTINUE
  590. 650 CONTINUE
  591. * on remet tout dans DHOOK
  592. DO 670 I=1,NSTRS
  593. DO 660 J=1,NSTRS
  594. DDHOOK(I,J)=GDMATR(I,J)
  595. 660 CONTINUE
  596. 670 CONTINUE
  597. ENDIF
  598. ENDIF
  599. *
  600. C Betocy
  601. ELSE IF (MAPL.EQ.54) THEN
  602. *
  603. * notice de l operateur MODE : "en 2D contraintes planes"
  604. * sinon on rend la matrice elastique
  605. IF (IFOU.EQ.-2) THEN
  606. CALL BETOKT(XMAT,TRAC,XVAR,XSTRS,NTRAT,NTRAC,DDHOOK,KERRE)
  607. * kerre toujours nul
  608. ENDIF
  609. *
  610. C Rotating Crack
  611. ELSE IF (MAPL.EQ.55) THEN
  612. *
  613. * notice de MATE : "la direction de fissuration change a chaque pas"
  614. * possible en 1D ? dans ROTKTA DDHOOK(4,4), en 1D on va jusque 3
  615. * sinon on rend la matrice elastique
  616. IF (IFOU.EQ.-2) THEN
  617. CALL ROTKTA(XMAT,XVAR,XSTRS,DDHOOK,KERRE)
  618. * kerre = 2 "impossible d inverser"
  619. IF (KERRE.EQ.2) IRET = 0
  620. ENDIF
  621. *
  622. C Plastique parfait, cinematique, isotrope
  623. ELSE IF ((MAPL.EQ.1).OR.(MAPL.EQ.4).OR.(MAPL.EQ.5)) THEN
  624. *
  625. IF (IFOU.EQ.-2) THEN
  626. * 2D contraintes planes
  627. XXXX = XSTRS(1)*XSTRS(1) + XSTRS(2)*XSTRS(2)
  628. YYYY = XSTRS(1)*XSTRS(2)
  629. ZZZZ = XSTRS(4)*XSTRS(4)
  630. AUX4=UNTIER/(UN-XNU*XNU)
  631. U=(DEUX*XSTRS(1)-XSTRS(2)+XNU*(DEUX*XSTRS(2)-
  632. & XSTRS(1)))*AUX4
  633. V=(DEUX*XSTRS(2)-XSTRS(1)+XNU*(DEUX*XSTRS(1)-
  634. & XSTRS(2)))*AUX4
  635. W= XSTRS(4)/(UN+XNU)
  636. A= CINQ*XXXX - HUIT*YYYY
  637. B= CINQ*YYYY - DEUX*XXXX
  638. C= DETIER*(XXXX-YYYY)+DEUX*ZZZZ
  639. D= (A+DEUX*XNU*B)*AUX4 + SIX*ZZZZ/(UN+XNU)
  640. XNUM=TROIS*YOU*(YOU-YOUTA)
  641. DENOM=DEUX*YOUTA*C + (YOU - YOUTA)*D
  642. BETA=XNUM/DENOM
  643. *
  644. DDHOOK(1,1)=DDHOOK(1,1)-BETA*U*U
  645. DDHOOK(2,1)=DDHOOK(2,1)-BETA*V*U
  646. DDHOOK(4,1)=DDHOOK(4,1)-BETA*W*U
  647. DDHOOK(1,2)=DDHOOK(1,2)-BETA*U*V
  648. DDHOOK(2,2)=DDHOOK(2,2)-BETA*V*V
  649. DDHOOK(4,2)=DDHOOK(4,2)-BETA*W*V
  650. DDHOOK(1,4)=DDHOOK(1,4)-BETA*U*W
  651. DDHOOK(2,4)=DDHOOK(2,4)-BETA*V*W
  652. DDHOOK(4,4)=DDHOOK(4,4)-BETA*W*W
  653. *
  654. ELSE IF ((IFOU.EQ.4).OR.(IFOU.EQ.8)) THEN
  655. * DYCZ,GYCZ
  656. XXXX = XSTRS(1)*XSTRS(1) + XSTRS(2)*XSTRS(2)
  657. YYYY = XSTRS(1)*XSTRS(2)
  658. AUX4=UNTIER/(UN-XNU*XNU)
  659. U=(DEUX*XSTRS(1)-XSTRS(2)+XNU*(DEUX*XSTRS(2)-
  660. & XSTRS(1)))*AUX4
  661. V=(DEUX*XSTRS(2)-XSTRS(1)+XNU*(DEUX*XSTRS(1)-
  662. & XSTRS(2)))*AUX4
  663. A= CINQ*XXXX - HUIT*YYYY
  664. B= CINQ*YYYY - DEUX*XXXX
  665. C= DETIER*(XXXX-YYYY)
  666. D= (A+DEUX*XNU*B)*AUX4
  667. XNUM=TROIS*YOU*(YOU-YOUTA)
  668. DENOM=DEUX*YOUTA*C + (YOU - YOUTA)*D
  669. BETA=XNUM/DENOM
  670. *
  671. DDHOOK(1,1)=DDHOOK(1,1)-BETA*U*U
  672. DDHOOK(2,1)=DDHOOK(2,1)-BETA*V*U
  673. DDHOOK(1,2)=DDHOOK(1,2)-BETA*U*V
  674. DDHOOK(2,2)=DDHOOK(2,2)-BETA*V*V
  675. *
  676. ELSE IF ((IFOU.EQ.5).OR.(IFOU.EQ.10).OR.(IFOU.EQ.13)
  677. & .OR.(IFOU.EQ.6)) THEN
  678. * CYDZ,CYGZ,AXCZ,CYCZ
  679. * comment adapter les formules 2D contraintes planes au 1D CY.. ?
  680. * pour le moment, on rend la matrice elastique
  681. *
  682. ELSE
  683. * tout sauf contraintes planes (1D aussi)
  684. * deviateur
  685. XMOY=(XSTRS(1)+XSTRS(2)+XSTRS(3))*UNTIER
  686. DO 201 IA=1,3
  687. DEVIAT(IA)=XSTRS(IA)- XMOY
  688. 201 CONTINUE
  689. IF (NSTRS.GE.4) THEN
  690. DO 301 IA=4,NSTRS
  691. DEVIAT(IA)=XSTRS(IA)
  692. 301 CONTINUE
  693. ENDIF
  694. *
  695. XNUM = XNEUF * YOU * ( YOU - YOUTA)
  696. XXXX = TROIS * ( YOU - YOUTA ) +
  697. & DEUX * ( UN + XNU ) * YOUTA
  698. DENOM= DEUX * SIGET2 * ( UN + XNU ) * XXXX
  699. BETA = XNUM / DENOM
  700. *
  701. DO 400 ICOMP1=1,NSTRS
  702. DO 400 ICOMP2=1,NSTRS
  703. DDHOOK(ICOMP1,ICOMP2)= DDHOOK(ICOMP1,ICOMP2) -
  704. & BETA*DEVIAT(ICOMP1)*DEVIAT(ICOMP2)
  705. 400 CONTINUE
  706. *
  707. ENDIF
  708. *
  709. C Drucker simple, Drucker Prager
  710. ELSE IF ((MAPL.EQ.3).OR.(MAPL.EQ.15)) THEN
  711. *
  712. * si SIGET = 0 (avec tr(S) dans le critere, c est possible)
  713. * on rend la matrice elastique (sinon on divise par SIGET=0)
  714. IF (SIGET.EQ.0.) RETURN
  715. CALL ZERO(DDF,NSTRS,1)
  716. CALL ZERO(DDG,NSTRS,1)
  717. *
  718. IF ((IFOU.EQ.-2).OR.(IFOU.EQ.4).OR.(IFOU.EQ.8)) THEN
  719. * 2D contraintes planes et 1D DYCZ GYCZ
  720. SIGM = (XSTRS(1)-XSTRS(2)*0.5D0)/SIGET
  721. DDF(1)=XMAT(1) + XMAT(2)*SIGM
  722. DDG(1)=XMAT(4) + XMAT(5)*SIGM
  723. SIGM = (XSTRS(2)-XSTRS(1)*0.5D0)/SIGET
  724. DDF(2)=XMAT(1) + XMAT(2)*SIGM
  725. DDG(2)=XMAT(4) + XMAT(5)*SIGM
  726. IF (NSTRS.GE.4) THEN
  727. DDF(4)=6.D0*XMAT(2)*XSTRS(4)/SIGET
  728. DDG(4)=6.D0*XMAT(5)*XSTRS(4)/SIGET
  729. ENDIF
  730. *
  731. ELSE IF ((IFOU.EQ.5).OR.(IFOU.EQ.10).OR.(IFOU.EQ.13)) THEN
  732. * 1D CYDZ CYGZ et AXCZ
  733. SIGM = (XSTRS(1)-XSTRS(3)*0.5D0)/SIGET
  734. DDF(1)=XMAT(1) + XMAT(2)*SIGM
  735. DDG(1)=XMAT(4) + XMAT(5)*SIGM
  736. SIGM = (XSTRS(3)-XSTRS(1)*0.5D0)/SIGET
  737. DDF(3)=XMAT(1) + XMAT(2)*SIGM
  738. DDG(3)=XMAT(4) + XMAT(5)*SIGM
  739. *
  740. ELSE IF (IFOU.EQ.6) THEN
  741. * 1D CYCZ
  742. SIGM = (XSTRS(1)-XSTRS(2)*0.5D0)/SIGET
  743. DDF(1)=XMAT(1) + XMAT(2)*SIGM
  744. DDG(1)=XMAT(4) + XMAT(5)*SIGM
  745. *
  746. ELSE
  747. SIGM = (XSTRS(1)-XSTRS(2)*0.5D0 -XSTRS(3)*0.5D0)/SIGET
  748. DDF(1)=XMAT(1) + XMAT(2)*SIGM
  749. DDG(1)=XMAT(4) + XMAT(5)*SIGM
  750. SIGM = (XSTRS(2)-XSTRS(1)*0.5D0 -XSTRS(3)*0.5D0)/SIGET
  751. DDF(2)=XMAT(1) + XMAT(2)*SIGM
  752. DDG(2)=XMAT(4) + XMAT(5)*SIGM
  753. SIGM = (XSTRS(3)-XSTRS(1)*0.5D0 -XSTRS(2)*0.5D0)/SIGET
  754. DDF(3)=XMAT(1) + XMAT(2)*SIGM
  755. DDG(3)=XMAT(4) + XMAT(5)*SIGM
  756. IF (NSTRS.GE.4) THEN
  757. DO 421 I=4,NSTRS
  758. DDF(I)=6.D0*XMAT(2)*XSTRS(I)/SIGET
  759. DDG(I)=6.D0*XMAT(5)*XSTRS(I)/SIGET
  760. 421 CONTINUE
  761. ENDIF
  762. *
  763. ENDIF
  764. *
  765. DO 422 I=1,NSTRS
  766. DHDF(I)=0.D0
  767. DHDG(I)=0.D0
  768. DO 422 J=1,NSTRS
  769. DHDF(I)=DHDF(I)+DDHOOK(I,J)*DDF(J)
  770. DHDG(I)=DHDG(I)+DDHOOK(I,J)*DDG(J)
  771. 422 CONTINUE
  772. *
  773. DENOM = XMAT(9)* SQRT(2.D0*XMAT(4)*XMAT(4)+
  774. & XMAT(5)*XMAT(5))
  775. DO 423 I=1,NSTRS
  776. DENOM = DENOM + DDF(I)*DHDG(I)
  777. 423 CONTINUE
  778. DO 320 ICOMP1=1,NSTRS
  779. DO 320 ICOMP2=1,NSTRS
  780. DDHOOK(ICOMP1,ICOMP2)=DDHOOK(ICOMP1,ICOMP2)
  781. & - DHDG(ICOMP1)*DHDF(ICOMP2)/DENOM
  782. 320 CONTINUE
  783. *
  784. ELSE
  785. * on rend la matrice elastique pour les autres modeles
  786. ENDIF
  787. *
  788. *-----------------------------------------------------------------------
  789. * MATERIAU UNIDIRECTIONNEL (MATE=4)
  790. * Acier unidirectionnel (MAPL=40)
  791. *-----------------------------------------------------------------------
  792. ELSE IF ((MATE.EQ.4).AND.(MAPL.EQ.40).AND.(IFOU.LE.0)) THEN
  793. *-----------------------------------------------------------------------
  794. DO 530 IC=1,3
  795. MPTVAL=IVAMAT
  796. MELVAL=IVAL(IC)
  797. IGMN=MIN(IGAU,VELCHE(/1))
  798. IBMN=MIN(IB ,VELCHE(/2))
  799. XMAT(IC)=VELCHE(IGMN,IBMN)
  800. 530 CONTINUE
  801. V1X=TXR(1,1)*XMAT(2)+TXR(1,2)*XMAT(3)
  802. V1Y=TXR(2,1)*XMAT(2)+TXR(2,2)*XMAT(3)
  803. XMAT(2)=V1X
  804. XMAT(3)=V1Y
  805. C Lecture Variables internes
  806. MPTVAL=IVARI
  807. MELVAL=IVAL(4)
  808. IGMN=MIN(IGAU,VELCHE(/1))
  809. IBMN=MIN(IB ,VELCHE(/2))
  810. XVAR(1)=VELCHE(IGMN,IBMN)
  811. *
  812. CALL UNIAKT(XMAT,XVAR,IFOU,DDHOOK,KERRE)
  813. * kerre toujours nul
  814. IRET = 1
  815. *
  816. *-----------------------------------------------------------------------
  817. * AUTRE MATERIAU :
  818. * UNIDIRECTIONNEL (sauf MAPL=40 et IFOU = -3,-2,-1,0)
  819. * ORTHOTROPE
  820. * ANISOTROPE
  821. *-----------------------------------------------------------------------
  822. ELSE
  823. *-----------------------------------------------------------------------
  824. IF (MATE.EQ.2) THEN
  825. CMATE = 'ORTHOTRO'
  826. ELSE IF (MATE.EQ.3) THEN
  827. CMATE = 'ANISOTRO'
  828. ELSE IF (MATE.EQ.4) THEN
  829. CMATE = 'UNIDIREC'
  830. ELSE
  831. C cet operateur ne fonctionne pas pour le materiau %m1:8
  832. IRET = -834
  833. RETURN
  834. ENDIF
  835. SEGINI,WRKMAT
  836. MPTVAL = IVAMAT
  837. SEGACT,MPTVAL
  838. DO 9004 IM=1,NMATT
  839. IF (IVAL(IM).NE.0) THEN
  840. MELVAL=IVAL(IM)
  841. IBMN=MIN(IB ,VELCHE(/2))
  842. IGMN=MIN(IGAU,VELCHE(/1))
  843. if (ibmn.gt.0.and.igmn.gt.0) then
  844. VALMAT(IM)= VELCHE(IGMN,IBMN)
  845. C* else
  846. C* VALMAT(IM)=0.D0
  847. endif
  848. C* ELSE
  849. C* VALMAT(IM)=0.D0
  850. ENDIF
  851. 9004 CONTINUE
  852. CALL DOHMAO(VALMAT,CMATE,IFOU,IDIM,TXR,XLOC,XGLOB,D1HO,
  853. & ROTH,DDHOOK,LHOOK,1,IRETOU)
  854. SEGSUP,WRKMAT
  855. IF (IRETOU.EQ.0) THEN
  856. IRET = -81
  857. ELSE IF (IRETOU.EQ.1) THEN
  858. IRET = 1
  859. ELSE
  860. IRET = 0
  861. ENDIF
  862. ENDIF
  863. *
  864. RETURN
  865. END
  866.  
  867.  
  868.  
  869.  

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