Télécharger syco12.eso

Retour à la liste

Numérotation des lignes :

  1. C SYCO12 SOURCE PV 17/12/08 21:17:49 9660
  2. C SYCO12 SOURCE FANDEUR 11/04/12 21:15:09 6938
  3.  
  4. SUBROUTINE SYCO12(wrk52,wrk53,wrk54,wrk2,wrk3,
  5. & IB,IGAU,NBPGAU,LTRAC,IFOURB,iecou,xecou)
  6.  
  7. C-----------------------------------------------------------------------
  8. C ECOULEMENT PLASTIQUE POUR UN POINT
  9. C ALGORITHME ORTIZ ET SIMO
  10. C
  11. C EN ENTREE :
  12. C
  13. C SIG0 CONTRAINTES AU DEBUT DU PAS
  14. C DEPST INCREMENT DE DEFORMATIONS TOTALES
  15. C ( THERMIQUE ENLEVEE )
  16. C VAR0 VARIABLES INTERNES DEDUT DU PAS
  17. C VAREX0 VARIABLES EXTERNES DEBUT DU PAS
  18. C VAREXF VARIABLES EXTERNES FIN DU PAS
  19. C XMAT COEFFICIENTS DU MATERIAU
  20. C PRECIS PRECISION POUR ITERATIONS INTERNES
  21. C WORK DES CARACTERISTIQUES
  22. C TRAC COURBE DE TRACTION
  23. C MFR1 INDICE DE FORMULATION
  24. C NSTRSS NOMBRE DE CONTRAINTES CA2000
  25. C INPLAS NUMERO DU MODELE DE PLASTICITE
  26. C DDAUX = MATRICE DE HOOKE ELASTIQUE
  27. C CMATE = NOM DU MATERIAU
  28. C VALMAT= TABLEAU DE CARACTERISTIQUES DU MATERIAU
  29. C VALCAR= TABLEAU DE CARACTERISTIQUES GEOMETRIQUES
  30. C N2EL = NBRE D ELEMENTS DANS SEGMENT DE HOOKE
  31. C N2PTEL= NBRE DE POINTS DANS SEGMENT DE HOOKE
  32. C IFOU = OPTION DE CALCUL
  33. C IB = NUMERO DE L ELEMENT COURANT
  34. C IGAU = NUMERO DU POINT COURANT
  35. C EPAIST= EPAISSEUR
  36. C NBPGAU= NBRE DE POINTS DE GAUSS
  37. C MELE = NUMERO DE L ELEMENT FINI
  38. C NPINT = NBRE DE POINTS D INTEGRATION
  39. C NBGMAT= NBRE DE POINTS DANS SEGMENT DE CARACTERISTIQUES
  40. C NELMAT= NBRE D ELEMENTS DANS SEGMENT DE CARACTERISTIQUES
  41. C SECT = SECTION
  42. C LHOOK = TAILLE DE LA MATRICE DE HOOKE
  43. C TXR,XLOC,XGLOB,D1HOOK,ROTHOO,DDHOMU,CRIGI = TABLEAUX UTILISES
  44. C UTILISES POUR LE CALCUL DE LA MATRICE DE HOOKE
  45. C
  46. C EN SORTIE :
  47. C
  48. C SIGF CONTRAINTES A LA FIN DU PAS
  49. C VARF VARIABLES INTERNES FIN DU PAS
  50. C DEFP INCR. DE DEFORMATIONS PLASTIQUES
  51. C KERRE CODE D'ERREUR
  52. C = 0 SI TOUT OK
  53. C = 99 SI FORMULATION NON DISPONIBLE
  54. C EN PLASTICITE
  55. C = 1 SI DEPS EST NEGATIF
  56. C = 2 SI NOMBRE MAX D'ITERATIONS INTERNES DEPASSE
  57. C
  58. C IFOUR : OPTION DE CALCUL
  59. C
  60. C IFOUR = -3 DEFORMATION PLANE GENERALISEE
  61. C IFOUR = -2 CONTRAINTES PLANES
  62. C IFOUR = -1 DEFORMATIONS PLANES
  63. C IFOUR = 0 AXISYMETRIQUE
  64. C IFOUR = 1 SERIE DE FOURIER
  65. C IFOUR = 2 TRIDIMENSIONNEL
  66. C
  67. C MFR1 : NUMERO DE LA FORMULATION ELEMENT FINI
  68. C
  69. C MFR1 = 1 MASSIF
  70. C MFR1 = 63 XFEM
  71. C
  72. C ISOTROPE
  73. C INPLAS : NUMERO DU MODELE DE VISCOPLASTICITE
  74. C
  75. C INPLAS = 153 SYCO1
  76. C INPLAS = 154 SYCO2
  77. C
  78. C-----------------------------------------------------------------------
  79. C CONVENTION DE REMPLISSAGE DES MEMOIRES
  80. C-----------------------------------------------------------------------
  81. C
  82. C XMAT(1) MODULE D'YOUNG
  83. C XMAT(2) COEFFICIENT DE POISSON
  84. C
  85. C TRAC(1 A 2*LTRAC) COURBE DE TRACTION
  86. C
  87. C MODELE ISOTROPE
  88. C VAR0(1) EPSE
  89. C VAR0(2) VP
  90. C
  91. C
  92. C-----------------------------------------------------------------------
  93. C 20/09/2017 : modif SG critere de convergence trop serre
  94. C TEST=PETI*APHI0 + utilisation CCREEL
  95. C idem ccoin0.eso
  96.  
  97. IMPLICIT INTEGER(I-N)
  98. IMPLICIT REAL*8(A-H,O-Z)
  99.  
  100. -INC CCOPTIO
  101. -INC CCREEL
  102. -INC DECHE
  103.  
  104. SEGMENT IECOU
  105. INTEGER icow1,icow2,icow3,icow4,icow5,icow6,icow7,
  106. 1 icow8,icow9,icow10,icow11,icow12,icow13,icow14,icow15,icow16,
  107. 2 icow17,icow18,icow19,icow20,icow21,icow22,icow23,icow24,
  108. 3 icow25,icow26,icow27,icow28,icow29,icow30,icow31,
  109. 4 icow32,icow33,NSTRSS,MFR1, NBGMAT,NELMAT,icow38,
  110. 5 icow39,icow40,icow41,icow42,icow43,icow44
  111. INTEGER icow45,icow46,icow47,icow48,icow49,icow50,
  112. . icow51,icow52,icow53,icow54,icow55,icow56
  113. . icow57,icow58
  114. ENDSEGMENT
  115. *
  116. SEGMENT WRK2
  117. REAL*8 TRAC(LTRAC)
  118. ENDSEGMENT
  119. *
  120. SEGMENT WRK3
  121. REAL*8 WORK(LW),WORK2(LW2)
  122. ENDSEGMENT
  123. *
  124. * Commun XECOU: sert de fourre-tout pour les initialisations
  125. * de réels
  126. *
  127. SEGMENT XECOU
  128. * COMMON/XECOU/DTOPTI,TSOM,TCAR,DTT,DT,TREFA,TEMP00
  129. REAL*8 DTOPTI,TSOM,TCAR,DTT,DT,TREFA,TEMP00
  130. ENDSEGMENT
  131. C COMMON/XECOU/DTOPTI,TSOM,TCAR,DTT,DT,TREFA,TEMP00
  132. C-----------------------------------------------------------------------
  133.  
  134. DIMENSION SIG(130),EPS(130)
  135. DIMENSION S(8),SX(8),DS(8),DSIG(8),SPHER(8),SPHER1(8),SPHER2(8)
  136. DIMENSION DSPHER1(8),DSPHER2(8),DEPSE(8),DEPSP(8),DDEPSE(8)
  137. DIMENSION F(8),W1(8),W2(8),SIGB(8),Z1(8),DIV(8),DDA(8,8)
  138. DIMENSION CRIGI(12)
  139.  
  140. NCONT=NSTRSS
  141. C
  142. C--------PAS DE TEMPS -----------------------
  143. TSYC9 = DT
  144. itracb=1
  145. 2222 continue
  146.  
  147. C-----------------------------------------------------------------------
  148. DO I=1,NSTRSS
  149. S(I)=0.D0
  150. SPHER(I)=0.D0
  151. SPHER1(I)=0.D0
  152. SPHER2(I)=0.D0
  153. ENDDO
  154. KERRE=0
  155. YUNG=XMAT(1)
  156. XNU=XMAT(2)
  157.  
  158. C---------CARACTERISTIQUES GEOMETRIQUES---------------------------------
  159. C
  160. C
  161. EP1=1.D0
  162. ALFAH=1.D0
  163.  
  164. IF(INPLAS.EQ.153) THEN
  165. C WRITE(6,*) 'INPLAS', inplas
  166. C WRITE(6,*) 'XMAT', xmat(1), xmat(2), xmat(3),xmat(4),xmat(5)
  167. C WRITE(6,*) 'XMAT', xmat(6), xmat(7), xmat(8)
  168. PSYC9=XMAT(6)
  169. DSYC9=XMAT(7)
  170. ENDIF
  171. IF(INPLAS.EQ.154) THEN
  172. C WRITE(6,*) 'INPLAS', inplas
  173. C WRITE(6,*) 'XMAT', xmat(1), xmat(2), xmat(3),xmat(4),xmat(5)
  174. C WRITE(6,*) 'XMAT', xmat(6), xmat(7), xmat(8),xmat(9),xmat(10)
  175. PSYC9=XMAT(6)
  176. ASYC9=XMAT(7)
  177. BSYC9=XMAT(8)
  178. CSYC9=XMAT(9)
  179. ENDIF
  180.  
  181. DO I=1,LTRAC
  182. SIG(I)=TRAC(2*I-1)
  183. EPS(I)=TRAC(2*I)
  184. ENDDO
  185.  
  186. EPST=VAR0(1)
  187. C
  188. C VP=VAR0(2)
  189. C if(VP.LE.1.d-4) VP = 1.d-4
  190. C VP0 = 1.D-4
  191. C IF(INPLAS.EQ.153) THEN
  192. C deno1 = 1.D0 + ( (VP0/DSYC9)**(1.D0/PSYC9) )
  193. C deno2 = 1.D0 + ( (VP/DSYC9)**(1.D0/PSYC9) )
  194. C else if(INPLAS.EQ.154) then
  195. C HSYC9 = ASYC9 + (BSYC9*EXP(-EPST/CSYC9))
  196. C deno1 = 1.D0 + (HSYC9 * ((VP0)**(1.D0/PSYC9)))
  197. C deno2 = 1.D0 + (HSYC9 * ((VP)**(1.D0/PSYC9)))
  198. C endif
  199. C
  200.  
  201. C-xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx-
  202. C-xxxx pour le materiau syco1 xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx-
  203. C DSYC9 parametre D de la loi Symonds Cowper
  204. C PSYC9 parametre P de la loi Symonds Cowper
  205. C VP valeur minimum de la vitesse de deformation equivalente
  206. IF(INPLAS.EQ.153) THEN
  207. VP0 = 1.D-4
  208. VP = VP0
  209. else if(INPLAS.EQ.154) then
  210. C VP valeur minimum de la vitesse de deformation equivalente
  211. VP0 = 1.D-4
  212. VP = VP0
  213. endif
  214. C-xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx-
  215. CALL CALSIG(DEPST,DDAUX,NSTRSS,CMATE,VALMAT,VALCAR,
  216. & N2EL,N2PTEL,MFR1,IFOURB,IB,IGAU,EPAIST,
  217. & NBPGAU,MELE,NPINT,NBGMAT,NELMAT,SECT,LHOOK,TXR,
  218. & XLOC,XGLOB,D1HOOK,ROTHOO,DDHOMU,CRIGI,DSIGT,IRTD)
  219.  
  220. IF(IRTD.NE.1) THEN
  221. KERRE=69
  222. GOTO 1000
  223. ENDIF
  224.  
  225. IF (ifourb.eq.-2) then
  226. * en cas de contraintes planes on repasse en 3D
  227. do iu=1,6
  228. do ju=1,6
  229. DDA(iu,ju)=0.d0
  230. enddo
  231. enddo
  232. cte_cp = xnu / (1.d0 - xnu)
  233. aux= YUNG / ((1.d0+xnu)*(1.d0-2.d0*xnu))
  234. aux1= aux * xnu
  235. aux2= aux * (1.d0-xnu)
  236. gege = Yung / (2.d0 *(1.d0 +xnu))
  237. DDA(1,1)=AUX2
  238. DDA(2,1)=AUX1
  239. DDA(1,2)=AUX1
  240. DDA(2,2)=AUX2
  241. DDA(3,3)=aux2
  242. DDA(1,3)=aux1
  243. DDA(2,3)=aux1
  244. DDA(3,1)=aux1
  245. DDA(3,2)=aux1
  246.  
  247. DDA(4,4)=gege
  248. DDA(5,5)=gege
  249. DDA(6,6)=GEGE
  250. ENDIF
  251. *
  252. DO I=1,NSTRSS
  253. S(I)=SIG0(I)+DSIGT(I)
  254. SIGB(I)=S(I)
  255. SX(I)=S(I)-SPHER(I)
  256. ENDDO
  257.  
  258.  
  259.  
  260. C-----------------------------------------------------------------------
  261. C CALCUL DE LA LIMITE ELASTIQUE SI
  262. C-----------------------------------------------------------------------
  263.  
  264. C
  265. CALL TRACTI(SI,EPST,SIG,EPS,LTRAC,2,IBI)
  266. c
  267. IF(IBI.EQ.1) THEN
  268. KERRE=75
  269. GOTO 1000
  270. ENDIF
  271.  
  272. C-----------------------------------------------------------------------
  273. C CALCUL DE LA CONTRAINTE EQUIVALENTE SEQ
  274. C-----------------------------------------------------------------------
  275. * attention en contraintes planes on se declare en 3D
  276. * rien besoin de faire dans vonmis0 car ifourb n'est pas utilisé
  277. * et les contraintes sont dimensionnées à 6
  278. * if(igau.eq.1) then
  279. * WRITE(6,*) ' en entrée sx ',sx(1),sx(2),sx(3)
  280. * endif
  281. C
  282. C SEQ=VONMIS0(SX,NSTRSS,MFR1,IFOURB,EP1,ALFAH)
  283. C calcul de la contrainte de Von Mises
  284. SOM1=0.D0
  285. DO I=4,NSTRSS
  286. SOM1=SOM1+SX(I)*SX(I)
  287. ENDDO
  288.  
  289. SEQ=SX(1)*SX(1)+SX(2)*SX(2)+SX(3)*SX(3)-SX(1)*SX(2)-SX(1)*SX(3)
  290. & -SX(2)*SX(3)+3.D0*SOM1
  291. SEQ=SQRT(ABS(SEQ))
  292. C
  293. C-----------------------------------------------------------------------
  294.  
  295. C LE CRITERE EST-IL VERIFIE
  296. C-----------------------------------------------------------------------
  297. C if(itracb.eq.1) then
  298. C WRITE(6,*)' en entrée si seq', si, seq
  299. C endif
  300. C
  301. c test de plasticite sur SI statique
  302. C SI = SI * deno2/deno1
  303. C
  304. PHI=SEQ-SI
  305. NITER=0
  306. PETI=1.1D0*0.5D0*PRECIS*SEQ
  307. CALL EPSPRE(SEQ,SI,PETI,ITRY)
  308. C
  309. IF((ITRY.EQ.1).OR.(SEQ.LE.SI)) THEN
  310. * rien a faire on n'a pas plastifié
  311. DO I=1,NCONT
  312. SIGF(I)=S(I)
  313. DEFP(I)=0.D0
  314. ENDDO
  315.  
  316. VARF(1)=VAR0(1)
  317. VARF(2)=VP0
  318. RETURN
  319. ENDIF
  320.  
  321. C-----------------------------------------------------------------------
  322. C ON A PLASTIFIE
  323. C-----------------------------------------------------------------------
  324. C
  325. cri0= si * 1.d-8
  326. PHI0=PHI
  327. SI0 = SI
  328. RR=0.D0
  329.  
  330. DO I=1,NCONT
  331. DSIG(I)=0.D0
  332. DEPSP(I)=0.D0
  333. DSPHER1(I)=0.D0
  334. DSPHER2(I)=0.D0
  335. ENDDO
  336.  
  337. C-----------------------------------------------------------------------
  338. C DEBUT DE LA BOUCLE D'ITERATIONS INTERNES
  339. C-----------------------------------------------------------------------
  340. sx1in=sx(1)
  341. sx2in=sx(2)
  342. sx3in= sx(3)
  343. s1in=s(1)
  344. s2in=s(2)
  345.  
  346. iderin=0
  347. 10 NITER=NITER+1
  348.  
  349. phi= seq - si
  350.  
  351.  
  352. C---------CALCUL DE W1=DF/D(SIGMA)--------------------------------------
  353.  
  354. C---------ELEMENTS MASSIFS----------------------------------------------
  355.  
  356.  
  357. F(1)=(2.D0*SX(1)-SX(2)-SX(3))/3.D0
  358. F(2)=(2.D0*SX(2)-SX(1)-SX(3))/3.D0
  359. F(3)=(2.D0*SX(3)-SX(1)-SX(2))/3.D0
  360. DO I=4,NSTRSS
  361. F(I)=SX(I)
  362. ENDDO
  363. DO I=1,3
  364. W1(I)=1.5D0*F(I)/SEQ
  365. Z1(I)=W1(I)
  366. ENDDO
  367. DO I=4,NCONT
  368. W1(I)=3.D0*F(I)/SEQ
  369. Z1(I)=1.5D0*F(I)/SEQ
  370. ENDDO
  371.  
  372.  
  373. if(ifourb.eq.-2) then
  374. DO I=1,NCONT
  375. r_z = 0.D0
  376. DO J=1,NCONT
  377. r_z=r_z+DDA(I,J)*W1(J)
  378. ENDDO
  379. W2(I)=r_z
  380. ENDDO
  381. else
  382. DO I=1,NCONT
  383. r_z = 0.D0
  384. DO J=1,NCONT
  385. r_z=r_z+DDAUX(I,J)*W1(J)
  386. ENDDO
  387. W2(I)=r_z
  388. ENDDO
  389. endif
  390. COEF=0.D0
  391. DO I=1,NCONT
  392. COEF=COEF+W1(I)*W2(I)
  393. ENDDO
  394.  
  395. C-----------------------------------------------------------------------
  396. C ECROUISSAGE ISOTROPE
  397. C-----------------------------------------------------------------------
  398.  
  399. CALL TRACTI(PENTE,EPST,SIG,EPS,LTRAC,1,IBI)
  400. * WRITE(6,*) ' pente epst', pente,epst
  401. IF(IBI.EQ.1) THEN
  402. KERRE=75
  403. GOTO 1000
  404. ENDIF
  405.  
  406. RP=PENTE
  407. C=0.D0
  408.  
  409. DENOM=COEF+C
  410. C-xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx-
  411. C---------CALCUL des termes en sup dus a la viscosité----------------
  412. C il s'agit de -qsi1 = -d(phi)/d(°p) = +d(SY)/d(°p)
  413. c (où °p = dp/dt cad la vitesse inelastiq equivalente)
  414. IF(INPLAS.EQ.153) THEN
  415. deno1 = 1.D0 + ( (VP0/DSYC9)**(1.D0/PSYC9) )
  416. deno2 = 1.D0 + ( (VP/DSYC9)**(1.D0/PSYC9) )
  417.  
  418. qsi1 = (SI/(PSYC9*DSYC9))*((VP/DSYC9)**( (1.D0/PSYC9)-1.D0))
  419. qsi1 = qsi1 / (deno2 * TSYC9)
  420. DENOM = DENOM + qsi1
  421. DENOM = DENOM + (RP * deno2/deno1)
  422. C
  423. ELSE IF (INPLAS.EQ.154) THEN
  424. HSYC9 = ASYC9 + (BSYC9*EXP(-EPST/CSYC9))
  425. deno1 = 1.D0 + (HSYC9 * ((VP0)**(1.D0/PSYC9)))
  426. deno2 = 1.D0 + (HSYC9 * ((VP)**(1.D0/PSYC9)))
  427.  
  428. qsi1 = (SI*HSYC9/PSYC9) * (VP**((1.D0/PSYC9)-1.D0))
  429. qsi1 = qsi1 / (deno2 * TSYC9)
  430. DENOM = DENOM + qsi1
  431. DENOM = DENOM + (RP * deno2/deno1)
  432. C
  433. qsi2 = (SI/deno2)*(VP**(1.D0/PSYC9))*((ASYC9-HSYC9)/CSYC9)
  434. DENOM = DENOM + qsi2
  435. C
  436. ENDIF
  437. C-xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx-
  438. C23456789012345678901234567890123456789012345678901234567890123456789012
  439. DELTA=PHI/DENOM
  440. DMU=C*DELTA/SEQ
  441. C if(itracb.eq.1) then
  442. C WRITE(6,*)'pente coef denom delta dmu niter si seq'
  443. C WRITE(6,*)pente ,coef , denom ,delta ,dmu,niter,si,seq
  444. C endif
  445.  
  446. DO I=1,NCONT
  447. DSIG(I)=-DELTA*W2(I)
  448. DSPHER1(I)=DMU*SX(I)
  449. ENDDO
  450.  
  451. * Cas des contraintes planes en massif
  452. if(ifourb.eq.-2) then
  453.  
  454. bb= abs(dsig(3)+ s(3) )
  455.  
  456. C if(itracb.eq.1) then
  457. C WRITE(6,*) ' iderin bb ' , iderin , bb
  458. C endif
  459.  
  460. r_z = dsig(3) * cte_cp
  461. sx(3)= sx3in - dsig(3)
  462. sx(1)= sx1in - r_z
  463. sx(2)= sx2in - r_z
  464. C
  465. C SEQ=VONMIS0(SX,NSTRSS,MFR1,IFOURB,EP1,ALFAH)
  466. C calcul de la contrainte de Von Mises
  467. SOM1=0.D0
  468. DO I=4,NSTRSS
  469. SOM1=SOM1+SX(I)*SX(I)
  470. ENDDO
  471.  
  472. SEQ=SX(1)*SX(1)+SX(2)*SX(2)+SX(3)*SX(3)-SX(1)*SX(2)-SX(1)*SX(3)
  473. & -SX(2)*SX(3)+3.D0*SOM1
  474. SEQ=SQRT(ABS(SEQ))
  475. C
  476. s(3)= - dsig(3)
  477. s(1)= s1in - r_z
  478. s(2)= s2in - r_z
  479.  
  480. cri0= si * 1.d-8
  481. if( bb.gt.cri0) then
  482.  
  483. if(iderin.eq.0) then
  484. niter=niter - 1
  485. endif
  486.  
  487. C if(itracb.eq.1) then
  488. C WRITE(6,* ) ' cri0 bb seq delta ',cri0,bb,seq , delta
  489. C endif
  490.  
  491. iderin = iderin + 1
  492.  
  493. ** WRITE(6,*) ' iderin ' , iderin
  494.  
  495. if(iderin.gt.50) then
  496. WRITE(6,*) 'on iderin superieur a 50'
  497. write(6,*) ' probleme dans iterations internes'
  498. KERRE=2
  499. GO TO 1000
  500. endif
  501.  
  502. C23456789012345678901234567890123456789012345678901234567890123456789012
  503. C WRITE(6,*)'cri0 bb niter iderin',cri0,bb,niter,iderin,dsig(3),s(3)
  504. go to 10
  505. endif
  506.  
  507. DMU=C*DELTA/SEQ
  508. DO I=1,NCONT
  509. DSPHER1(I)=DMU*SX(I)
  510. ENDDO
  511. endif
  512.  
  513. iderin=0
  514. C
  515. DR=RP*DELTA
  516. RR=RR+DR
  517. EPST=EPST+DELTA
  518.  
  519. C mise à jour des contraintes
  520. DO I=1,NCONT
  521. S(I)=S(I)+DSIG(I)
  522. SPHER1(I)=SPHER1(I)+DSPHER1(I)
  523. SPHER2(I)=SPHER2(I)+DSPHER2(I)
  524. SPHER(I)=SPHER1(I)+SPHER2(I)
  525. SX(I)=S(I)-SPHER(I)
  526. ENDDO
  527. if(ifourb.eq.-2) then
  528. s(3)=0.d0
  529. endif
  530. C if(igau.eq.1) then
  531. C WRITE(6,*) ' niter',niter
  532. C WRITE(6,*) ' dsig' , dsig(1),dsig(2),dsig(3),dsig(4)
  533. C WRITE(6,*) ' s' , s(1),s(2),s(3),s(4)
  534. C WRITE(6,*) ' sx' , sx(1),sx(2),sx(3),sx(4)
  535. C WRITE(6,*) ' spher' , spher(1),spher(2),spher(3),spher(4)
  536. C endif
  537.  
  538. C SEQ=VONMIS0(SX,NSTRSS,MFR1,IFOURB,EP1,ALFAH)
  539. C calcul de la contrainte de Von Mises
  540. SOM1=0.D0
  541. DO I=4,NSTRSS
  542. SOM1=SOM1+SX(I)*SX(I)
  543. ENDDO
  544.  
  545. SEQ=SX(1)*SX(1)+SX(2)*SX(2)+SX(3)*SX(3)-SX(1)*SX(2)-SX(1)*SX(3)
  546. & -SX(2)*SX(3)+3.D0*SOM1
  547. SEQ=SQRT(ABS(SEQ))
  548. C
  549. C---------CONTRAINTES PLANES--------------------------------------------
  550.  
  551. IF(IFOURB.EQ.-2) THEN
  552.  
  553. F(1)=(2.D0*SX(1)-SX(2)-SX(3))/3.D0
  554. F(2)=(2.D0*SX(2)-SX(1)-SX(3))/3.D0
  555. F(3)=(2.D0*SX(3)-SX(1)-SX(2))/3.D0
  556. DO I=4,NSTRSS
  557. F(I)=SX(I)
  558. ENDDO
  559. DO I=1,3
  560. W1(I)=1.5D0*F(I)/SEQ
  561. ENDDO
  562. DO I=4,NSTRSS
  563. W1(I)=3.D0*F(I)/SEQ
  564. ENDDO
  565.  
  566.  
  567. DO I=1,NSTRSS
  568. DEPSP(I)=DEPSP(I)+DELTA*W1(I)
  569. ENDDO
  570. ENDIF
  571.  
  572. C-----------------------------------------------------------------------
  573. C TEST
  574. C CALCUL DE LA NOUVELLE VALEUR DE PHI
  575. C-----------------------------------------------------------------------
  576. CALL TRACTI(SI,EPST,SIG,EPS,LTRAC,2,IBI)
  577. C******************************************************
  578. c-xxxx calcul de la valeur (i+1) de H, de la vitesse de def VP,
  579. c et de SI en prenant en compte la viscosité
  580. IF(INPLAS.EQ.153) THEN
  581. C materiau SYCO1
  582. deno1 = 1.D0 + ( (VP0/DSYC9)**(1.D0/PSYC9) )
  583. VP = MAX( (VP+(DELTA/TSYC9)) , VP0 )
  584. C deno2 = 1.D0 + ( (VP/DSYC9)**(1.D0/PSYC9) )
  585. SI = SI * (1.D0 + ( (VP/DSYC9)**(1.D0/PSYC9) )) / deno1
  586. C
  587. else if (INPLAS.EQ.154) then
  588. C materiau SYCO2
  589. HSYC9 = ASYC9 + (BSYC9*EXP(-EPST/CSYC9))
  590. deno1 = 1.D0 + (HSYC9 * ((VP0)**(1.D0/PSYC9)))
  591. VP = MAX( (VP+(DELTA/TSYC9)) , VP0 )
  592. C deno2 = 1.D0 + (HSYC9 * ((VP)**(1.D0/PSYC9)))
  593. SI = SI * (1.D0 + (HSYC9 * (VP**(1.D0/PSYC9)))) / deno1
  594. endif
  595. C******************************************************
  596.  
  597.  
  598. PHI=SEQ-SI
  599.  
  600. PETI=1.D-7
  601. APHI=ABS(PHI)
  602. APHI0=ABS(PHI0)
  603. TEST=max(PETI*APHI0,XZPREC*100.D0*SEQ)
  604. *sg TEST=PETI*APHI0
  605.  
  606.  
  607. IF(NITER.GT.200) THEN
  608. WRITE(6,*) 'on a bien planté ici,car APHI=',APHI
  609. WRITE(6,*) 'SI0,PHI0,RR,SI,PHI',SI0,PHI0,RR,SI,PHI
  610. TEST2=max(1.D-6*APHI0,XZPREC*100.D0*SEQ)
  611. *sg TEST2 = 1.D-6*APHI0
  612. if(APHI.gt.TEST2) then
  613. KERRE=2
  614. GO TO 1000
  615. else
  616. TEST=TEST2
  617. endif
  618. ELSEIF(NITER.eq.101) THEN
  619. WRITE(6,*) 'on aurait du planté ici,car APHI=',APHI
  620. WRITE(6,*) 'SI0,PHI0,RR,SI,PHI',SI0,PHI0,RR,SI,PHI
  621. WRITE(6,*) 'mais on va continuer 50pas de plus'
  622. ENDIF
  623.  
  624. * if(igau.eq.1) then
  625. * WRITE(6,*)'test conver niter phi seq si testphi0',niter,
  626. * $ phi,seq,si,test,aphi0
  627. * WRITE(6,*) ' rr si0' , rr, si0
  628. * endif
  629.  
  630. IF(APHI.LE.TEST) THEN
  631.  
  632.  
  633. C---------TOUTES FORMULATIONS SAUF CONTRAINTES PLANES-------------------
  634.  
  635. IF(IFOURB.NE.-2) THEN
  636. DO I=1,NCONT
  637. DS(I)=S(I)-SIGB(I)
  638. ENDDO
  639. C CALL EPSIG0(DS,DDEPSE,MFR1,IFOURB,YUNG,XNU,WORK,NSTRSS)
  640. C copie d'un bout de EPSIG0
  641. DO I=1,NSTRSS
  642. DDEPSE(I)=0.D0
  643. ENDDO
  644. AUX1=1.D0/YUNG
  645. AUX2=2.D0*(1.D0+XNU)
  646. AUX3=AUX1*AUX2
  647. DDEPSE(1)=(DS(1)-XNU*(DS(2)+DS(3)))*AUX1
  648. DDEPSE(2)=(DS(2)-XNU*(DS(1)+DS(3)))*AUX1
  649. DDEPSE(3)=(DS(3)-XNU*(DS(1)+DS(2)))*AUX1
  650. DO I=4,NSTRSS
  651. DDEPSE(I)=AUX3*DS(I)
  652. ENDDO
  653. C fin copie d'un bout de EPSIG0
  654. DO I=1,NCONT
  655. DEPSE(I)=DEPST(I)+DDEPSE(I)
  656. DEPSP(I)=DEPST(I)-DEPSE(I)
  657. ENDDO
  658. ENDIF
  659.  
  660. DO I=1,NSTRSS
  661. SIGF(I)=S(I)
  662. DEFP(I)=DEPSP(I)
  663. ENDDO
  664.  
  665. VARF(1)=EPST
  666. VARF(2)=VP
  667.  
  668. KERRE=0
  669. RETURN
  670.  
  671. ELSE
  672. sx1in=sx(1)
  673. sx2in=sx(2)
  674. s1in=s(1)
  675. s2in=s(2)
  676.  
  677. GOTO 10
  678. ENDIF
  679.  
  680. 1000 CONTINUE
  681. RETURN
  682. END
  683.  
  684.  
  685.  
  686.  
  687.  
  688.  
  689.  
  690.  
  691.  
  692.  
  693.  
  694.  

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