Télécharger deriv1.eso

Retour à la liste

Numérotation des lignes :

  1. C DERIV1 SOURCE AM 09/12/15 21:15:15 6587
  2. SUBROUTINE DERIV1(TAU,SIG,EPSV,VAR,SIGD,EPSVD,VARD,EPSTHD,
  3. &DSPT,XMAT,NSTRS,NVARI,NCOMAT,DD,DDV,DDINV,DDINVp,YOUNG,NYOUNG,
  4. &XNU,NXNU,MFR,XCAR,ICARA,IFOURB,INDIC,TI,TPOINT,ITHHER,VARF,
  5. & IB,IGAU,KERRE)
  6.  
  7. C-----------------------------------------------------------------------
  8. C Objet: Cette subroutine calcule les derivees des variables internes
  9. C pour un materiau viscoplastique a endommagement et ecrouissage
  10. C isotrope en regime anisotherme
  11. C-----------------------------------------------------------------------
  12. C
  13. C-----------------------------------------------------------------------
  14. C Entree: TAU pas de temps
  15. C TI temperature au debut du pas
  16. C TPOINT derivee de la temperature
  17. C SIG(NSTRS) tenseur des contraintes dans le repère général
  18. C EPSV(NSTRS) tenseur des deformations dans le repère général
  19. C VAR(NVARI) tableau contenant les variables internes
  20. C pilotant les equations
  21. C VAR(1)=p ; VAR(2)=D11 ; VAR(3)=D22; VAR(4)=D33;
  22. C VAR(5)=D12; VAR(6)=D13; VAR(7)=D23; VAR(8)=DMAX
  23. C EPSTHD vitesse de dilatation thermique au debut du pas
  24. C DSPT(NSTRS,NSTRS) vitesse de deformation totale
  25. C XMAT(NCOMAT) tableau des parametres scalaires du materiau
  26. C a une temperature TI donnee
  27. C XMAT(1)=YOUNG ;XMAT(2)=XNU ;XMAT(3)=N
  28. C XMAT(4)=M ;XMAT(5)=KK; XMAT(6)=ALF2;
  29. C XMAT(7)=BET2 ;XMAT(8)=r; XMAT(9)=A ;
  30. C XMAT(10)=q ;XMAT(11)=ALPHAT
  31. C XMAT(12)=RHO; XMAT(13)=SIGY
  32. C MFR indice de la formulation mecanique(seulement massif
  33. C ou coque pour les materiaux endommageables
  34. C ICARA nombre de caracteristiques geometriques des elements
  35. C finis
  36. C XCAR(ICARA) tableau des caracteristiques geometriques des
  37. C elements finis
  38. C IFOURB= -2 EN CONTR. PLANES
  39. C -1 EN DEFORM. PLANES
  40. C 0 EN AXISYMETRIE
  41. C 1 EN SERIE DE FOURIER
  42. C 2 EN TRIDIM
  43. C INDIC=0, 1 OU -1 pour plasticite avec endommagement
  44. C =2 OU -2 pour viscoplasticite avec endommagement
  45. C ITHHER = 0 pas de chargement thermique et materiau constant
  46. C = 1 chargement thermique et materiau constant
  47. C = 2 chargement thermique et materiau(T)
  48. C------------------------------------------------------------------------
  49. C
  50. C------------------------------------------------------------------------
  51. C Sortie: EPSVD(NSTRS) derivee du tenseur des deformations
  52. C VARD(NVARI) tableau contenant les derivees des variables
  53. C internes
  54. C VARF(NVARI) état des variables internes à la fin du sous programme
  55. C SIGD(NSTRS) derivee des contraintes
  56. C DD(NSTRS,NSTRS) matrice de Hooke au debut du pas
  57. C DDV(NSTRS,NSTRS) derivee de DD
  58. C DDINV(NSTRS,NSTRS) inverse de DD
  59. C------------------------------------------------------------------------
  60. C
  61. C
  62. IMPLICIT INTEGER(I-N)
  63. IMPLICIT REAL*8(A-H,O-Z)
  64. DIMENSION SIG(*),EPSV(*),VAR(*),XCAR(*),YOUNG(*),XNU(*)
  65. DIMENSION SIGD(*),EPSVD(*),VARD(*),DDINV(NSTRS,NSTRS)
  66. DIMENSION XMAT(*),DSPT(*),EPSTHD(*),VARF(NVARI)
  67. DIMENSION DDV(NSTRS,NSTRS),DDINVp(NSTRS,NSTRS),DD(NSTRS,NSTRS)
  68. DIMENSION EPSED(6),XX(6),Di(3),VaP(3),AA(3),BB(6)
  69. DIMENSION CO(3,3),COM(3,3),DCOT(3,3),COE(3,3),DCOE(3,3)
  70. DIMENSION AAA(3,3),BBB(3,3),CCC(3,3),UN2(3,3)
  71. DIMENSION FTH(3,3),ROT(3,3),ROTT(3,3),VcP(3,3),S(3,3)
  72. DIMENSION D1(3,3),D2(3,3),DV1(3,3),DI1(3,3),DI2(3,3)
  73. PARAMETER (UN=1.D0, UNM=-1.D0)
  74. LOGICAL endonul
  75. C
  76. C
  77. endonul = .FALSE.
  78. IF ((VAR(2).EQ.0).AND.(VAR(3).EQ.0).AND.(VAR(4).EQ.0).AND.
  79. & (VAR(5).EQ.0).AND.(VAR(6).EQ.0).AND.(VAR(7).EQ.0).AND.
  80. & (VAR(8).EQ.0)) THEN
  81. endonul = .TRUE.
  82. ENDIF
  83.  
  84. C
  85. C-----------------------------------------------------------------
  86. C
  87. C EXPRESSIONS DES TENSEURS UTILISES DANS LE MODELE
  88. C
  89. C-----------------------------------------------------------------
  90. C
  91. C----------------------------------------------
  92. C TENSEUR DES CONTRAINTES (dans repere general)
  93. C----------------------------------------------
  94. C
  95. CALL ZERO(CO,3,3)
  96. IF (mfr.eq.1) THEN
  97. DO I=1,3
  98. CO(I,I) = SIG(I)
  99. ENDDO
  100. CO(1,2) = SIG(4)
  101. IF(NSTRS.NE.6) GOTO 50
  102. CO(1,3) = SIG(5)
  103. CO(2,3) = SIG(6)
  104. 50 CONTINUE
  105. CO(2,1) = CO(1,2)
  106. CO(3,1) = CO(1,3)
  107. CO(3,2) = CO(2,3)
  108. ENDIF
  109. C
  110. C------------------------------------
  111. C TENSEUR D ENDOMMAGEMENT ET
  112. C TENSEUR D ENDOMMAGEMENT DIAGONALISE
  113. C------------------------------------
  114. C
  115. XD11=VAR(2)
  116. XD22=VAR(3)
  117. XD33=VAR(4)
  118. XD12=VAR(5)
  119. XD13=VAR(6)
  120. XD23=VAR(7)
  121. CALL ZERO(D1,3,3)
  122. D1(1,1) = XD11
  123. D1(2,2) = XD22
  124. D1(3,3) = XD33
  125. D1(1,2) = XD12
  126. D1(1,3) = XD13
  127. D1(2,3) = XD23
  128. D1(2,1) = D1(1,2)
  129. D1(3,1) = D1(1,3)
  130. D1(3,2) = D1(2,3)
  131.  
  132. *
  133. * AM 14/12/09 TEST SELON NSTRS
  134. *
  135. IF(IFOURB.LE.0) THEN
  136. JDIM = 2
  137. CALL JACOB4(D1,JDIM,Di,ROT)
  138. ELSE
  139. JDIM = 3
  140. CALL JACOB4(D1,JDIM,Di,ROT)
  141. ENDIF
  142. *
  143. CALL ZERO(D2,3,3)
  144. D2(1,1) = Di(1)
  145. D2(2,2) = Di(2)
  146. D2(3,3) = Di(3)
  147. w1 = 1.D0 - Di(1)
  148. w2 = 1.D0 - Di(2)
  149. w3 = 1.D0 - Di(3)
  150. IF (D2(1,1).GE.1.D0) w1=1.D-6
  151. IF (D2(2,2).GE.1.D0) w2=1.D-6
  152. IF (D2(3,3).GE.1.D0) w3=1.D-6
  153.  
  154. ADMAX = DMAX1(Di(1),Di(2),Di(3))
  155. C
  156. C------------------------
  157. C TENSEUR UNITE D ORDRE 2
  158. C------------------------
  159. C
  160. CALL ZERO(UN2,3,3)
  161. DO I=1,3
  162. UN2(I,I) = 1.D0
  163. ENDDO
  164. C
  165. C-------------------------------------------------------
  166. C TENSEUR DES CONTRAINTES MOTRICES (dans repere general)
  167. C-------------------------------------------------------
  168. C
  169. EQ=XMAT(10)
  170. DO I=1,3
  171. DO J=1,3
  172. ROTT(I,J)=ROT(I,J)
  173. ENDDO
  174. ENDDO
  175. CALL INVALM (ROTT, 3 , 3 , KER, 1.D-12)
  176. C
  177. IF (endonul) THEN
  178. DO I=1,3
  179. DO J=1,3
  180. COM(I,J) = CO(I,J)
  181. ENDDO
  182. ENDDO
  183. ELSE
  184. C --- Matrice de [(I-D)/det(I-D)]**(q/2) dans le repère principal de D ---
  185. CALL ZERO(DI1,3,3)
  186. DO I=1,3
  187. XNB=(1.D0-Di(I))/(w1*w2*w3)
  188. IF (XNB.GE.0.D0) THEN
  189. DI1(I,I)=(XNB)**(EQ/2.D0)
  190. ELSE
  191. DI1(I,I)=0.D0
  192. WRITE(6,*) 'erreur1, IB = ', IB, ' Igau=',igau
  193. KERRE = 99
  194. GO TO 100
  195. c Do j = 1,7
  196. c write(6,*) 'var(',j,')=', var(j)
  197. c ENDDO
  198. ENDIF
  199. ENDDO
  200. C --- passage dans le repère de général ---
  201. CALL MULMAT(AAA,ROT,DI1,3,3,3)
  202. CALL MULMAT(DI2,AAA,ROTT,3,3,3)
  203.  
  204. CALL MULMAT(AAA,DI2,CO,3,3,3)
  205. CALL MULMAT(COM,AAA,DI2,3,3,3)
  206. ENDIF
  207. C
  208. C-----------------------------------
  209. C TENSEUR DES CONTRAINTES EFFECTIVES
  210. C-----------------------------------
  211. C
  212. IF (endonul) THEN
  213. DO I=1,3
  214. DO J=1,3
  215. COE(I,J) = CO(I,J)
  216. ENDDO
  217. ENDDO
  218. ELSE
  219. C --- Matrice de [I-D]**(-1/2) dans le repère principal de D ---
  220. CALL ZERO(DI1,3,3)
  221. CALL ZERO(DI2,3,3)
  222. DO I=1,3
  223. XNB=1.D0-Di(I)
  224. IF (XNB.GT.0.D0) THEN
  225. DI1(I,I)=SQRT(1.D0/(XNB))
  226. ELSE
  227. DI1(I,I)=0.D0
  228. WRITE(6,*) 'erreur2 IB=', IB, ' IGAU= ', IGAU
  229. KERRE=99
  230. GO TO 100
  231. ENDIF
  232. ENDDO
  233. C --- passage dans le repère de général ---
  234. CALL MULMAT(AAA,ROT,DI1,3,3,3)
  235. CALL MULMAT(DI2,AAA,ROTT,3,3,3)
  236. c DO i=1,3
  237. c DO j=1,3
  238. c WRITE(6,*) 'DI2(',i,',',j,')=',DI2(i,j)
  239. c enddo
  240. c enddo
  241.  
  242. CALL MULMAT(AAA,DI2,CO,3,3,3)
  243. CALL MULMAT(COE,AAA,DI2,3,3,3)
  244. ENDIF
  245. C
  246. C----------------------------------
  247. C DEVIATEUR DE CONTRAINTE EFFECTIVE
  248. C----------------------------------
  249. C
  250. CALL ZERO(BBB,3,3)
  251. DO I=1,3
  252. BBB(I,I)=-(COE(1,1)+COE(2,2)+COE(3,3))/(3.D0)
  253. ENDDO
  254. CALL PLMATR(COE,UN,BBB,UN,3,3,DCOE)
  255. C --- Calcul de (I-D)**(-1/2).DCOE.(I-D)**(-1/2) ---
  256. IF (endonul) THEN
  257. DO I=1,3
  258. DO J=1,3
  259. DCOT(I,J) = DCOE(I,J)
  260. ENDDO
  261. ENDDO
  262. ELSE
  263. CALL MULMAT(AAA,DI2,DCOE,3,3,3)
  264. CALL MULMAT(DCOT,AAA,DI2,3,3,3)
  265. ENDIF
  266. C
  267. C----------------------------------------------------------------------
  268. C
  269. C CALCUL DE J1,J2
  270. C J1: CONTRAINTE MOYENNE
  271. C J2: CONTRAINTE DE VON MISES
  272. C
  273. C----------------------------------------------------------------------
  274. C
  275. AJ1 = (SIG(1)+SIG(2)+SIG(3))/(3.D0)
  276. DO I=1,3
  277. XX(I) = SIG(I)-AJ1
  278. ENDDO
  279. DO I=4,NSTRS
  280. XX(I) = SIG(I)
  281. ENDDO
  282. AJ2 = PROCON(XX,XX,NSTRS)
  283. AJ2 = SQRT(1.5D0*AJ2)
  284. C
  285. C----------------------
  286. C FORCE THERMODYNAMIQUE
  287. C----------------------
  288. C
  289. C --- Etat de triaxialité de la contrainte ---
  290. IF ((SIG(1).EQ.0).AND.(SIG(2).EQ.0).AND.(SIG(3).EQ.0)) THEN
  291. tri=1.D0
  292. ELSE
  293. tri = (AJ2)/(SIG(1)+SIG(2)+SIG(3))
  294. ENDIF
  295. c write(6,*) 'AJ2= ',AJ2
  296. c write(6,*) 'SIG(1)= ',SIG(1)
  297. c write(6,*) 'SIG(2)= ',SIG(2)
  298. c write(6,*) 'SIG(3)= ',SIG(3)
  299.  
  300. c write(6,*) 'tri= ',tri
  301.  
  302. C --- Valeurs propres et vecteurs propres du tenseur des contraintes motrices ---
  303. CALL JACOB3(COM,3,VaP,VcP)
  304. ER=XMAT(8)
  305. EA=XMAT(9)
  306. DO I=1,3
  307. AA(I) = ((tri*(VaP(I))/(EA))+ABS(tri*(VaP(I))/(EA)))/(2.D0)
  308. AA(I) = (AA(I))**(ER)
  309. ENDDO
  310. c DO I=1,3
  311. c DO J=1,3
  312. c write(6,*) 'COM(',i,',',j,')= ',COM(I,J)
  313. c ENDDO
  314. c ENDDO
  315. C DO J=1,3
  316. C write(6,*) 'VaP(',j,')= ',VaP(J)
  317. C ENDDO
  318. c DO J=1,3
  319. c write(6,*) 'AA(',j,')= ',AA(J)
  320. c ENDDO
  321. c DO I=1,3
  322. c DO J=1,3
  323. c write(6,*) 'VcP(',i,',',j,')= ',VcP(I,J)
  324. c ENDDO
  325. c ENDDO
  326. C
  327. C -- Dans le repère général ---
  328. CALL ZERO(FTH,3,3)
  329. FTH(1,1) = AA(1)*(VcP(1,1))**2+AA(2)*(VcP(1,2))**2
  330. FTH(1,1) = FTH(1,1)+AA(3)*(VcP(1,3))**2
  331. FTH(2,2) = AA(1)*(VcP(2,1))**2+AA(2)*(VcP(2,2))**2
  332. FTH(2,2) = FTH(2,2)+AA(3)*(VcP(2,3))**2
  333. FTH(3,3) = AA(1)*(VcP(3,1))**2+AA(2)*(VcP(3,2))**2
  334. FTH(3,3) = FTH(3,3)+AA(3)*(VcP(3,3))**2
  335. FTH(1,2) = AA(1)*VcP(1,1)*VcP(2,1)+AA(2)*VcP(1,2)*VcP(2,2)
  336. FTH(1,2) = FTH(1,2)+AA(3)*VcP(1,3)*VcP(2,3)
  337. FTH(1,3) = AA(1)*VcP(1,1)*VcP(3,1)+AA(2)*VcP(1,2)*VcP(3,2)
  338. FTH(1,3) = FTH(1,3)+AA(3)*VcP(1,3)*VcP(3,3)
  339. FTH(2,3) = AA(1)*VcP(2,1)*VcP(3,1)+AA(2)*VcP(2,2)*VcP(3,2)
  340. FTH(2,3) = FTH(1,3)+AA(3)*VcP(2,3)*VcP(3,3)
  341. FTH(3,1) = FTH(1,3)
  342. FTH(2,1) = FTH(1,2)
  343. FTH(3,2) = FTH(2,3)
  344. c DO I=1,3
  345. c DO J=1,3
  346. c write(6,*) 'FTH(',i,',',i,')= ',FTH(I,i)
  347. c ENDDO
  348. c ENDDO
  349. C
  350. C---------------------------------------------------------------------
  351. C
  352. C CALCUL DE LA VITESSE DE DEFORMATION PLASTIQUE CUMULEE VARD(1)
  353. C
  354. C---------------------------------------------------------------------
  355. C
  356. C --- Paramètres matériau ---
  357. EM=XMAT(4)
  358. EN=XMAT(3)
  359. EK=XMAT(5)
  360. ALF2=XMAT(6)
  361. C --- Calcul de la contrainte équivalente effective ---
  362. AJ3 = (COE(1,1)+COE(2,2)+COE(3,3))/(3.D0)
  363. DO I=1,3
  364. XX(I) = COE(I,I)-AJ3
  365. ENDDO
  366. XX(4) = COE(1,2)
  367. XX(5) = COE(1,3)
  368. XX(6) = COE(2,3)
  369. AJ4 = PROCON(XX,XX,NSTRS)
  370. AJ4 = SQRT(1.5D0*AJ4)
  371. C
  372. C --- Calcul de la norme euclidienne ---
  373. DO I=1,3
  374. DO J=1,3
  375. AAA(I,J)=(DCOT(I,J)*3.D0)/(2.D0*AJ4)
  376. ENDDO
  377. ENDDO
  378. CALL PLMATR(UN2,UN,D1,UNM,3,3,BBB)
  379. CALL INVALM (BBB, 3 , 3 , KER, 1.D-12)
  380. DO I=1,3
  381. DO J=1,3
  382. BBB(I,J)=BBB(I,J)*ALF2/3.D0
  383. ENDDO
  384. ENDDO
  385. CALL ZERO(CCC,3,3)
  386. CALL PLMATR(AAA,UN,BBB,UN,3,3,CCC)
  387. c DO I=1,3
  388. c DO j=1,3
  389. c write(6,*) 'CCC(',i,',',j,')= ', CCC(i,j)
  390. c ENDDO
  391. c ENDDO
  392.  
  393. val = 0.D0
  394. DO I=1,3
  395. DO J=1,3
  396. val = val+(CCC(I,J))**(2.D0)
  397. ENDDO
  398. ENDDO
  399. xnor = SQRT(val)
  400. C WRITE(6,*) 'xnor=',xnor
  401. C --- Calcul du premier terme du produit ---
  402. IF ( VAR(1).EQ.0.D0 ) VAR(1)=VAR(1)+1.D-12
  403. ter = (AJ4+(ALF2*AJ3))/(EK*(VAR(1)**(1.D0/EM)))
  404. ter = ter + ABS((AJ4+(ALF2*AJ3))/(EK*(VAR(1)**(1.D0/EM))))
  405. ter = ter/(2.D0)
  406. ter = ter**(EN)
  407. ter = ter*SQRT(2.D0/3.D0)
  408. c WRITE(6,*) 'ter=',ter
  409. C --- Calcul de l evolution de p ---
  410. IF ((SIG(1).EQ.0).AND.(SIG(2).EQ.0).AND.(SIG(3).EQ.0)) THEN
  411. VARD(1) = 0.D0
  412. ELSE
  413. VARD(1) = ter*xnor
  414. ENDIF
  415. C write(6,*) 'vard(1)=', vard(1)
  416. C
  417. C---------------------------------------------------------------------
  418. C
  419. C CALCUL DE LA VARIATION DE L ENDOMMAGEMENT VARD(I) (I=2..7)
  420. C
  421. C---------------------------------------------------------------------
  422.  
  423. BET2=XMAT(7)
  424. C --- Dans le repère principal de D ---
  425. CALL ZERO(DV1,3,3)
  426. DV1(1,1) = (BET2*FTH(1,1))+FTH(2,2)+FTH(3,3)
  427. DV1(2,2) = (BET2*FTH(2,2))+FTH(1,1)+FTH(3,3)
  428. DV1(3,3) = (BET2*FTH(3,3))+FTH(2,2)+FTH(1,1)
  429. DV1(1,2) = (BET2-1.D0)*FTH(1,2)
  430. DV1(1,3) = (BET2-1.D0)*FTH(1,3)
  431. DV1(2,3) = (BET2-1.D0)*FTH(2,3)
  432. DV1(2,1) = DV1(1,2)
  433. DV1(3,1) = DV1(1,3)
  434. DV1(3,2) = DV1(2,3)
  435. C
  436. C --- Expressions des VARD(I) (I=2..7) ---
  437. DO I=2,4
  438. VARD(I) = DV1(I-1,I-1)
  439. ENDDO
  440. VARD(5) = DV1(1,2)
  441. VARD(6) = DV1(1,3)
  442. VARD(7) = DV1(2,3)
  443. c DO I=2,4
  444. c write(6,*) 'vard(2)= ', vard(2)
  445. c ENDDO
  446. C
  447. C---------------------------------------------------------------------
  448. C
  449. C CALCUL DE LA VARIATION DE DEFORMATION PLASTIQUE
  450. C
  451. C---------------------------------------------------------------------
  452.  
  453. IF ((SIG(1).EQ.0).AND.(SIG(2).EQ.0).AND.(SIG(3).EQ.0)) THEN
  454. DO I=1,NSTRS
  455. EPSVD(I) = 0.D0
  456. ENDDO
  457. ELSE
  458. DO I=1,3
  459. EPSVD(I) = CCC(I,I)*VARD(1)*SQRT(3.D0/2.D0)/xnor
  460. ENDDO
  461. EPSVD(4) = CCC(1,2)*VARD(1)*SQRT(3.D0/2.D0)/xnor
  462. *
  463. * AM CORRECTION LE 9/12/09
  464. *
  465. **** IF (IFOUR.GT.0) THEN
  466. IF (IFOURB.GT.0) THEN
  467. EPSVD(5) = CCC(1,3)*VARD(1)*SQRT(3.D0/2.D0)/xnor
  468. EPSVD(6) = CCC(2,3)*VARD(1)*SQRT(3.D0/2.D0)/xnor
  469. ENDIF
  470. ENDIF
  471.  
  472.  
  473. c DO I=1,3
  474. c write(6,*) 'EPSVD(',i,')= ', EPSVD(i)
  475. c ENDDO
  476. C
  477. C----------------------------------------------------------------------
  478. C
  479. C CALCUL DE LA VARIATION DES CONTRAINTES
  480. C
  481. C----------------------------------------------------------------------
  482.  
  483. CALL DERTRA(NYOUNG,YOUNG,TI,YUNG,YUNGV,TINF,TSUP)
  484. CALL DERTRA(NXNU,XNU,TI,ENU,ENUV,TINF,TSUP)
  485. C
  486. C --- Calcul de l inverse de la matrice de Hooke au debut du pas ---
  487. CALL ELAST4(2,IFOURB,VAR,NVARI,XMAT,NCOMAT,YUNGV,ENUV,
  488. &XCAR,ICARA,MFR,NSTRS,DDINV,DDV,KERRE,INDIC,ITHHER)
  489. IF(KERRE.NE.0) GO TO 100
  490. C
  491. C --- Calcul de la matrice de Hooke et de sa derivee au debut du pas ---
  492. CALL ELAST4(1,IFOURB,VAR,NVARI,XMAT,NCOMAT,YUNGV,ENUV,
  493. &XCAR,ICARA,MFR,NSTRS,DD,DDV,KERRE,INDIC,ITHHER)
  494. IF(KERRE.NE.0) GO TO 100
  495. C
  496. C --- Détermination de la variation de la matrive de Hooke inversée ---
  497. CALL HOOKP(VAR,VARD,NVARI,XMAT,NCOMAT,MFR,NSTRS,DDINVp)
  498. C
  499. C --- Calcul de la variation de contrainte ---
  500. CALL ZDANUL(SIGD,NSTRS)
  501. DO I=1,3
  502. XX(I) = CO(I,I)
  503. ENDDO
  504. XX(4) = CO(1,2)
  505. XX(5) = CO(1,3)
  506. XX(6) = CO(2,3)
  507.  
  508. CALL ZDANUL(EPSED,NSTRS)
  509. DO I=1,NSTRS
  510. EPSED(I)=DSPT(I)-EPSVD(I)-EPSTHD(I)
  511. C write(6,*) 'DSPT(',I,') =', DSPT(I)
  512. c write(6,*) 'EPSVD(',I,') =', EPSVD(I)
  513. c write(6,*) 'EPSTHD(',I,') =', EPSTHD(I)
  514. C write(6,*) 'EPSED(',I,') =', EPSED(I)
  515. ENDDO
  516.  
  517. CALL MULMAT(BB,DDINVp,XX,NSTRS,1,NSTRS)
  518.  
  519. DO I=1,NSTRS
  520. XX(I)=EPSED(I)-BB(I)
  521. ENDDO
  522.  
  523. CALL MULMAT(SIGD,DD,XX,NSTRS,1,NSTRS)
  524. c DO I=1,3
  525. c IF (iGAU.EQ.1) THEN
  526. c write(6,*) 'SIGD(1) =', SIGD(1)
  527. c ENDIF
  528. c ENDDO
  529. C
  530. C----------------------------------------------------------------------
  531. C
  532. C MISE A JOUR DES VARIABLES INTERNES
  533. C
  534. C----------------------------------------------------------------------
  535.  
  536. VARF(8)=ADMAX
  537.  
  538. 100 CONTINUE
  539. RETURN
  540. END
  541.  
  542.  
  543.  
  544.  
  545.  

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