Télécharger codjf2.eso

Retour à la liste

Numérotation des lignes :

codjf2
  1. C CODJF2 SOURCE PV 17/12/08 21:16:15 9660
  2. C MODJF2 SOURCE AM 00/12/13 21:41:10 4045
  3. SUBROUTINE CODJF2(WRK52,WRK53,EPSP1,EPSP2,YLIM1,YLIM2,Z1,Z2,Y1,Y2,
  4. 1 JSIGN,SIG1,SIG2,CONV,DT1)
  5. C
  6. C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  7. C % %
  8. C % INTEGRATION DU MODJFD2 UNILATERAL AVEC FERMETURE LINEAIRE %
  9. C % DE FISSURES EN CONTRAINTES PLANES %
  10. C % (MODCLB NON ASSOCIE) %
  11. C % A L"AIDE DE LA METHODE DE DICHOTOMIE %
  12. C % AJOUT VICOENDOMMAGEMENT J.F. DUBE %
  13. C % %
  14. C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  15. C
  16. C
  17. C
  18. IMPLICIT INTEGER(I-N)
  19. -INC DECHE
  20. C VARIABLES D"ENTREES
  21. C
  22. C
  23. REAL*8 EPSP1,EPSP2
  24. C
  25. C VARIABLES DE SORTIE
  26. C
  27. C
  28. REAL*8 SIG1,SIG2,Z1,Z2,Y1,Y2,DT1
  29. INTEGER JSIGN(3)
  30. LOGICAL CONV
  31. C
  32. C
  33. C PARAMETRES D"EVOLUTION
  34. C
  35. C
  36. REAL*8 E,NU,A1,A2,B1,B2,Y01,Y02,SIGF1
  37. REAL*8 BETA1,BETA2,YLIM1,YLIM2
  38. C
  39. C
  40. C VARIABLES INTERNES
  41. C
  42. C
  43. REAL*8 Y001,Y002,DY01,DY02,hdp1,hdp2
  44. REAL*8 DELTASIG
  45. INTEGER NDIM
  46. PARAMETER (NDIM=6)
  47. REAL*8 JACOB(NDIM,NDIM),J15,J16,J25,J26
  48. REAL*8 DELAFON(NDIM)
  49. REAL*8 X(NDIM),DX(NDIM),X0(NDIM),XX,XIG1,XIG2
  50. REAL*8 ERREUR,TEST
  51. REAL*8 H1,H2,HT,HT2,trf
  52. REAL*8 ALFA1,ALFA2
  53. REAL*8 SIGY
  54. REAL*8 UN,DEUX
  55. REAL*8 ZMAX
  56. INTEGER ITER ,I,DIMPB
  57. LOGICAL LDOM1,LDOM2,evol1,evol2
  58. REAL*8 M1,M2,ALPHA1,ALPHA2
  59. REAL*8 x3,x4,x5,x6,dx5,dx6
  60. INTEGER IMAX
  61. REAL*8 pas1,erreurd1,xd1,xg1,pas2,erreurd2,xd2,xg2
  62. REAL*8 erreur_max,fctd1,fctd2
  63. logical passage1,passage2
  64. UN=1.D0
  65. DEUX=2.D0
  66. ZMAX=1.d+04
  67. E = XMAT(1)
  68. NU = XMAT(2)
  69. Y001 = XMAT(5)
  70. Y002 = XMAT(6)
  71. A1 = XMAT(7)
  72. A2 = XMAT(8)
  73. B1 = XMAT(9)
  74. B2 = XMAT(10)
  75. BETA1 = XMAT(11)
  76. BETA2 = XMAT(12)
  77. SIGF1 = XMAT(13)
  78. M1 = XMAT(14)
  79. M2 = XMAT(15)
  80. ALPHA1= XMAT(16)
  81. ALPHA2= XMAT(17)
  82. SIGY =SQRT(2.D0*E*Y001)
  83. * DELTASIG=SIGY*1.D-4
  84. DELTASIG=0.d0
  85. * print*,'E,NU,Y01,Y02,A1,A2,B1,B2,BETA1,BETA2'
  86. * print*,E,NU,Y01,Y02,A1,A2,B1,B2,BETA1,BETA2
  87. TEST=1.D-12
  88. ERREUR=1.D0
  89. ITER=0
  90. X(1)=SIG1
  91. X(2)=SIG2
  92. X(3)=Y1
  93. X(4)=Y2
  94. X(5)=Z1
  95. X(6)=Z2
  96. DO I=1, NDIM
  97. X0(I)=X(I)
  98. END DO
  99. x5=Z1
  100. x6=z2
  101. C
  102. C
  103. C
  104. JACOB(1,1) = -1.d0
  105. JACOB(1,2) = 0.D0
  106. JACOB(1,3) = 0.D0
  107. JACOB(1,4) = 0.D0
  108. JACOB(1,5) = 0.D0
  109. JACOB(1,6) = 0.D0
  110. JACOB(2,1) = 0.D0
  111. JACOB(2,2) = -1.d0
  112. JACOB(2,3) = 0.D0
  113. JACOB(2,4) = 0.D0
  114. JACOB(2,5) = 0.D0
  115. JACOB(2,6) = 0.D0
  116. JACOB(3,1) = 0.D0
  117. JACOB(3,2) = 0.D0
  118. JACOB(3,3) = -1.d0
  119. JACOB(3,4) = 0.D0
  120. JACOB(3,5) = 0.D0
  121. JACOB(3,6) = 0.D0
  122. JACOB(4,1) = 0.D0
  123. JACOB(4,2) = 0.D0
  124. JACOB(4,3) = 0.D0
  125. JACOB(4,4) = -1.d0
  126. JACOB(4,5) = 0.D0
  127. JACOB(4,6) = 0.D0
  128. JACOB(5,1) = 0.D0
  129. JACOB(5,2) = 0.D0
  130. JACOB(5,3) = 0.D0
  131. JACOB(5,4) = 0.D0
  132. JACOB(5,5) = -1.d0
  133. JACOB(5,6) = 0.D0
  134. JACOB(6,1) = 0.D0
  135. JACOB(6,2) = 0.D0
  136. JACOB(6,3) = 0.D0
  137. JACOB(6,4) = 0.D0
  138. JACOB(6,5) = 0.D0
  139. JACOB(6,6) = -1.d0
  140. Erreur_max=1.0d-10
  141. pas1=1.0
  142. Xd1=z1
  143. Xg1=z1
  144. imax=0
  145. Erreurd1=1.0
  146. passage1=.true.
  147. Do while ( (Erreurd1 .GT. Erreur_max) .and. (pas1 .ge. 1.d-6) )
  148. x(5)=xd1+pas1
  149.  
  150. c print*,'*** Z1 ***'
  151. c print*,x(5)
  152.  
  153. ** Calcul de Y01 modifie par l'effet de vitesse
  154. c dp1=(DIM(x(5),z1)/(x(5)*z1*dt1))
  155. c y01=y001+m1*dp1**alpha1
  156.  
  157. pas2=1.0
  158. Xd2=z2
  159. Xg2=x(6)
  160. Erreurd2=1.0
  161. passage2=.true.
  162. X(4)=Y2
  163. Do while ((Erreurd2 .GT. Erreur_max) .and. (pas2 .ge. 1.d-6))
  164. if (imax .ge. 10000) then
  165. print*,'pb dans jfd2'
  166. conv=.false.
  167. return
  168. end if
  169. x(6)=xd2+pas2
  170.  
  171. c print*,'*** Z2 ***'
  172. c print*,x(6)
  173.  
  174. ** Calcul de Y02 modifie par l'effet de vitesse
  175. c dp2=(DIM(x(6),z2)/(x(6)*z2*dt1))
  176. c y02=y002+m2*dp2**alpha2
  177. ** Calcul de sigma et Y1 et Y2
  178. ERREUR=1.D0
  179. ITER=0
  180. X(1)=Sig1
  181. X(2)=Sig2
  182. DO WHILE (ERREUR .GT. TEST .AND. ITER .LE. 100)
  183. ITER=ITER+1
  184. IF (X(1) .GE. DELTASIG) THEN
  185. H1=UN
  186. ELSE IF (X(1) .LE. -DELTASIG) THEN
  187. H1=0.D0
  188. ELSE
  189. H1=(X(1)+DELTASIG)/(DEUX*DELTASIG)
  190. END IF
  191. IF (X(2) .GE. DELTASIG) THEN
  192. H2=UN
  193. ELSE IF (X(2) .LE. -DELTASIG) THEN
  194. H2=0.D0
  195. ELSE
  196. H2=(X(2)+DELTASIG)/(DEUX*DELTASIG)
  197. END IF
  198. c
  199. c smooth transition on ht ht2
  200. c
  201. trf = x(1) + x(2)
  202. IF(trf.GE.deltasig)THEN
  203. HT=UN
  204. HT2=0.D0
  205. else if(trf.ge.(-deltasig)) then
  206. ht = (trf+deltasig)/DEUX/deltasig
  207. ht2 = UN-ht
  208. else if(trf.ge.(-sigf+deltasig)) then
  209. ht = 0.D+00
  210. ht2 = UN
  211. else if(trf.ge.(-sigf-deltasig)) then
  212. ht = 0.D+00
  213. c ht2 =UN- (trf-(sigf-deltasig))/(DEUX*sigf)
  214. ht2 = (trf+sigf+deltasig)/DEUX/deltasig
  215. else
  216. ht = 0.D+00
  217. ht2 = 0.D+00
  218. endif
  219. c print*,'DTjfd',dt1
  220. ** Calcul de Y01 modifie par l'effet de vitesse
  221. dy01=(DIM(x(5),z1)/(x(5)*z1*dt1))
  222. hdp1=1.0
  223. if (dy01 .gt. 0.0) hdp1=0.0
  224. y01=y001+m1*dy01**alpha1
  225. ** Calcul de Y02 modifie par l'effet de vitesse
  226. dy02=(DIM(x(6),z2)/(x(6)*z2*dt1))
  227. hdp2=1.0
  228. if (dy02 .gt. 0.0) hdp2=0.0
  229. y02=y002+m2*dy02**alpha2
  230. c print*,'dans clbeco iter=',iter,' y01=',y01,'y02=',y02
  231.  
  232. ldom1=.false.
  233. evol1=.false.
  234. ldom2=.false.
  235. evol2=.false.
  236. c print*,'y1',x(3),'ylim1',ylim1
  237. c print*,'z1',x(5),'x5',x5
  238. c print*,'ldom1',ldom1,'ldom2',ldom2
  239. C =====> PROBLEME ELASTIQUE
  240. C
  241. C
  242. C CALCUL DU JACOBIEN
  243. C
  244. C
  245. JACOB(3,1) = H1*X(1)*X(5)*X(5)/E
  246. C
  247. JACOB(3,2) = H2*X(2)*X(5)*X(5)/E
  248. C
  249. JACOB(4,1) = (UN-H1)*X(1)*X(6)*X(6)/E
  250. C
  251. JACOB(4,2) = (UN-H2)*X(2)*X(6)*X(6)/E
  252. C
  253. C CALCUL DU SECOND MEMBRE
  254. C
  255. XIG1 =((((X(5)-UN)*BETA1*HT2+(X(5)-UN)*BETA1*HT-E*EPSP2+(X(6)-UN)
  256. 1 *BETA2)*NU+BETA1*(X(6)*(H2+X(5)*(UN-H2)-UN)+X(5)**2*H2-X(5)*H2)
  257. 2 *HT2+BETA1*(X(6)*(H2+X(5)*(UN-H2)-UN)+X(5)**2*H2-X(5)*H2)*HT+
  258. 3 BETA2*(X(6)*(X(5)*H2+H2-UN)-X(5)*H2+X(6)**2*(UN-H2))+E*EPSP1*
  259. 4 (X(6)*(H2-UN)-X(5)*H2))*SIGF1+(BETA1*E*(X(5)*EPSP2-EPSP2)+
  260. 5 (UN-X(5))*BETA1*E*EPSP1)*HT2)/((NU**2+H1*(X(5)*X(6)*(2*H2-UN)-
  261. 6 X(5)**2*H2+X(6)**2*(UN-H2))-X(5)*X(6)*H2+X(6)**2*(H2-UN))*SIGF1+
  262. 7 (DEUX-DEUX*X(5))*BETA1*HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*
  263. 8 (-H2+X(5)*(H2-DEUX)+DEUX))+((X(5)-UN)*X(6)-X(5)**2+X(5))
  264. 9 *BETA1*H1)*HT2)
  265. C
  266. XIG2 =((((X(5)-UN)*BETA1*HT2+(X(5)-UN)*BETA1*HT-E*EPSP1+(X(6)-UN)
  267. 1 *BETA2)*NU+(((UN-X(5))*X(6)+X(5)**2-X(5))*BETA1*H1+(X(5)-UN)*
  268. 2 X(6)*BETA1)*HT2+(((UN-X(5))*X(6)+X(5)**2-X(5))*BETA1*H1+
  269. 3 (X(5)-UN)*X(6)*BETA1)*HT+(E*(X(6)*EPSP2-X(5)*EPSP2)+(-X(6)**2+
  270. 4 (X(5)+UN)*X(6)-X(5))*BETA2)*H1-X(6)*E*EPSP2+(X(6)**2-X(6))*
  271. 5 BETA2)*SIGF1+(BETA1*E*(EPSP2-X(5)*EPSP2)+(X(5)-UN)*BETA1*E*EPSP1)
  272. 6 *HT2)/((NU**2+H1*(X(5)*X(6)*(DEUX*H2-UN)-X(5)**2*H2+X(6)**2*
  273. 7 (UN-H2))-X(5)*X(6)*H2+X(6)**2*(H2-UN))*SIGF1+(DEUX-DEUX*X(5))*
  274. 8 BETA1*HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*
  275. 9 (H2-DEUX)+DEUX))+((X(5)-UN)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)
  276.  
  277. DELAFON(1) =X(1)-XIG1
  278. C
  279. DELAFON(2) =X(2)-XIG2
  280. C
  281. XX= (H1*X(1)*X(1)+H2*X(2)*X(2))*X(5)*X(5)/(DEUX*E)
  282. DELAFON(3) = X(3)-XX
  283. C
  284. XX=((UN-H1)*X(1)*X(1)+(UN-H2)*X(2)*X(2) )
  285. 1 *X(6)*X(6)/(DEUX*E)
  286. DELAFON(4) = X(4)-XX
  287. C
  288. IF (LDOM1) THEN
  289. C
  290. C =====> ENDOMMAGEMENT EN TRACTION
  291. C
  292. J15 = (((BETA1*HT2+BETA1*HT)*NU+BETA1*(DEUX*X(5)*H2-H2+X(6)*(UN-H2
  293. & ))*HT2+BETA1*(DEUX*X(5)*H2-H2+X(6)*(UN-H2))*HT+BETA2*(X(6)*H2-H
  294. & 2)-E*EPSP1*H2)*SIGF1+(BETA1*E*EPSP2-BETA1*E*EPSP1)*HT2)/((NU**2+
  295. & H1*(X(5)*X(6)*(DEUX*H2-UN)-X(5)**2*H2+X(6)**2*(UN-H2))-X(5)*X(6
  296. & )*H2+X(6)**2*(H2-UN))*SIGF1+(DEUX-DEUX*X(5))*BETA1*HT2*NU+(BETA1
  297. & *(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-DEUX)+DEUX))+((X(5)-UN
  298. & )*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)-((H1*(X(6)*(DEUX*H2-UN)-DEU
  299. & X*X(5)*H2)-X(6)*H2)*SIGF1-DEUX*BETA1*HT2*NU+(BETA1*(-DEUX*X(5)*H
  300. & 2+H2+X(6)*(H2-DEUX))+(X(6)-DEUX*X(5)+UN)*BETA1*H1)*HT2)*((((X(5
  301. & )-UN)*BETA1*HT2+(X(5)-UN)*BETA1*HT-E*EPSP2+(X(6)-UN)*BETA2)*NU+
  302. & BETA1*(X(6)*(H2+X(5)*(UN-H2)-UN)+X(5)**2*H2-X(5)*H2)*HT2+BETA1*
  303. & (X(6)*(H2+X(5)*(UN-H2)-UN)+X(5)**2*H2-X(5)*H2)*HT+BETA2*(X(6)*(
  304. & X(5)*H2+H2-UN)-X(5)*H2+X(6)**2*(UN-H2))+E*EPSP1*(X(6)*(H2-UN)-X
  305. & (5)*H2))*SIGF1+(BETA1*E*(X(5)*EPSP2-EPSP2)+(UN-X(5))*BETA1*E*EPS
  306. & P1)*HT2)/((NU**2+H1*(X(5)*X(6)*(DEUX*H2-UN)-X(5)**2*H2+X(6)**2*
  307. & (UN-H2))-X(5)*X(6)*H2+X(6)**2*(H2-UN))*SIGF1+(DEUX-DEUX*X(5))*BE
  308. & TA1*HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-DEUX)
  309. & +DEUX))+((X(5)-UN)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)**2
  310. C
  311. C
  312. J25 = (((BETA1*HT2+BETA1*HT)*NU+((-X(6)+DEUX*X(5)-UN)*BETA1*H1+X(6
  313. & )*BETA1)*HT2+((-X(6)+DEUX*X(5)-UN)*BETA1*H1+X(6)*BETA1)*HT+((X(
  314. & 6)-UN)*BETA2-E*EPSP2)*H1)*SIGF1+(BETA1*E*EPSP1-BETA1*E*EPSP2)*HT
  315. & 2)/((NU**2+H1*(X(5)*X(6)*(DEUX*H2-UN)-X(5)**2*H2+X(6)**2*(UN-H2
  316. & ))-X(5)*X(6)*H2+X(6)**2*(H2-UN))*SIGF1+(DEUX-DEUX*X(5))*BETA1*HT
  317. & 2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-DEUX)+DEUX)
  318. & )+((X(5)-UN)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)-((H1*(X(6)*(DEUX
  319. & *H2-UN)-DEUX*X(5)*H2)-X(6)*H2)*SIGF1-DEUX*BETA1*HT2*NU+(BETA1*(-
  320. & DEUX*X(5)*H2+H2+X(6)*(H2-DEUX))+(X(6)-DEUX*X(5)+UN)*BETA1*H1)*H
  321. & T2)*((((X(5)-UN)*BETA1*HT2+(X(5)-UN)*BETA1*HT-E*EPSP1+(X(6)-UN)
  322. & *BETA2)*NU+(((UN-X(5))*X(6)+X(5)**2-X(5))*BETA1*H1+(X(5)-UN)*X(
  323. & 6)*BETA1)*HT2+(((UN-X(5))*X(6)+X(5)**2-X(5))*BETA1*H1+(X(5)-UN)
  324. & *X(6)*BETA1)*HT+(E*(X(6)*EPSP2-X(5)*EPSP2)+(-X(6)**2+(X(5)+UN)*
  325. & X(6)-X(5))*BETA2)*H1-X(6)*E*EPSP2+(X(6)**2-X(6))*BETA2)*SIGF1+(B
  326. & ETA1*E*(EPSP2-X(5)*EPSP2)+(X(5)-UN)*BETA1*E*EPSP1)*HT2)/((NU**2
  327. & +H1*(X(5)*X(6)*(DEUX*H2-UN)-X(5)**2*H2+X(6)**2*(UN-H2))-X(5)*X(
  328. & 6)*H2+X(6)**2*(H2-UN))*SIGF1+(DEUX-DEUX*X(5))*BETA1*HT2*NU+(BETA
  329. & 1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-DEUX)+DEUX))+((X(5)-U
  330. & N)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)**2
  331. C
  332. JACOB(1,5) = J15
  333. JACOB(2,5) = J25
  334. JACOB(3,5) = (H1*X(1)*X(1)+H2*X(2)*X(2))*X(5)/E
  335. JACOB(4,5) = 0.D0
  336. JACOB(5,5) = -1.D0
  337. & -A1*B1*(A1*(X(3)-Y01))**(B1-UN)
  338. & *alpha1*m1*dy01**alpha1/(dy01+hdp1)/(x(5)**2*dt1)
  339. JACOB(5,1) = 0.D0
  340. JACOB(5,2) = 0.D0
  341. JACOB(5,3) = A1*B1*(A1*(X(3)-Y01))**(B1-UN)
  342. JACOB(5,4) = 0.D0
  343. DELAFON(5) = X(5)-((A1*(X(3)-Y01))**B1+UN)
  344.  
  345. IF (LDOM2) THEN
  346. C
  347. C =====> ENDOMMAGEMENT EN COMPRESSION
  348. C
  349. J16 = (BETA2*NU+BETA1*(H2+X(5)*(1-H2)-1)*HT2+BETA1*(H2+X(5)*(1-H2)
  350. 1 -1)*HT+BETA2*(X(5)*H2+H2+2*X(6)*(1-H2)-1)+E*EPSP1*(H2-1))*SIGF1/
  351. 2 ((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X(5)**2*H2+X(6)**2*(1-H2))-X(5)*
  352. 3 X(6)*H2+X(6)**2*(H2-1))*SIGF1+(2-2*X(5))*BETA1*HT2*NU+(BETA1*(-X
  353. 4 (5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-2)+2))+((X(5)-1)*X(6)-X(5)
  354. 5 **2+X(5))*BETA1*H1)*HT2)-((H1*(X(5)*(2*H2-1)+2*X(6)*(1-H2))-X(5
  355. 6 )*H2+2*X(6)*(H2-1))*SIGF1+(BETA1*(-H2+X(5)*(H2-2)+2)+(X(5)-1)*BE
  356. 7 TA1*H1)*HT2)*((((X(5)-1)*BETA1*HT2+(X(5)-1)*BETA1*HT-E*EPSP2+(X
  357. 8 (6)-1)*BETA2)*NU+BETA1*(X(6)*(H2+X(5)*(1-H2)-1)+X(5)**2*H2-X(5)
  358. 9 *H2)*HT2+BETA1*(X(6)*(H2+X(5)*(1-H2)-1)+X(5)**2*H2-X(5)*H2)*HT+
  359. : BETA2*(X(6)*(X(5)*H2+H2-1)-X(5)*H2+X(6)**2*(1-H2))+E*EPSP1*(X(6
  360. ; )*(H2-1)-X(5)*H2))*SIGF1+(BETA1*E*(X(5)*EPSP2-EPSP2)+(1-X(5))*BE
  361. < TA1*E*EPSP1)*HT2)/((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X(5)**2*H2+X(6
  362. = )**2*(1-H2))-X(5)*X(6)*H2+X(6)**2*(H2-1))*SIGF1+(2-2*X(5))*BETA1
  363. > *HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-2)+2))+(
  364. ? (X(5)-1)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)**2
  365. C
  366. J26 = (BETA2*NU+((1-X(5))*BETA1*H1+(X(5)-1)*BETA1)*HT2+((1-X(5))*B
  367. 1 ETA1*H1+(X(5)-1)*BETA1)*HT+(E*EPSP2+(-2*X(6)+X(5)+1)*BETA2)*H1-
  368. 2 E*EPSP2+(2*X(6)-1)*BETA2)*SIGF1/((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X
  369. 3 (5)**2*H2+X(6)**2*(1-H2))-X(5)*X(6)*H2+X(6)**2*(H2-1))*SIGF1+(2-
  370. 4 2*X(5))*BETA1*HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)
  371. 5 *(H2-2)+2))+((X(5)-1)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)-((H1*(X
  372. 6 (5)*(2*H2-1)+2*X(6)*(1-H2))-X(5)*H2+2*X(6)*(H2-1))*SIGF1+(BETA1*
  373. 7 (-H2+X(5)*(H2-2)+2)+(X(5)-1)*BETA1*H1)*HT2)*((((X(5)-1)*BETA1*H
  374. 8 T2+(X(5)-1)*BETA1*HT-E*EPSP1+(X(6)-1)*BETA2)*NU+(((1-X(5))*X(6)
  375. 9 +X(5)**2-X(5))*BETA1*H1+(X(5)-1)*X(6)*BETA1)*HT2+(((1-X(5))*X(6
  376. : )+X(5)**2-X(5))*BETA1*H1+(X(5)-1)*X(6)*BETA1)*HT+(E*(X(6)*EPSP2
  377. ; -X(5)*EPSP2)+(-X(6)**2+(X(5)+1)*X(6)-X(5))*BETA2)*H1-X(6)*E*EPS
  378. < P2+(X(6)**2-X(6))*BETA2)*SIGF1+(BETA1*E*(EPSP2-X(5)*EPSP2)+(X(5)
  379. = -1)*BETA1*E*EPSP1)*HT2)/((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X(5)**2*
  380. > H2+X(6)**2*(1-H2))-X(5)*X(6)*H2+X(6)**2*(H2-1))*SIGF1+(2-2*X(5))
  381. ? *BETA1*HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-2)
  382. @ +2))+((X(5)-1)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)**2
  383. C
  384. C
  385. JACOB(1,6) = J16
  386. JACOB(2,6) = J26
  387. JACOB(3,6) = 0.D0
  388. JACOB(4,6) = ((1.D0-H1)*X(1)*X(1)+(1.D0-H2)
  389. & *X(2)*X(2))*X(6)/E
  390. JACOB(5,6) = 0.D0
  391. JACOB(6,6) = -1.D0
  392. & -A2*B2*(A2*(X(4)-Y02))**(B2-UN)
  393. & *alpha2*m2*dy02**alpha2/(dy02+hdp2)/(x(6)**2*dt1)
  394. JACOB(6,1) = 0.D0
  395. JACOB(6,2) = 0.D0
  396. JACOB(6,3) = 0.D0
  397. * JACOB(6,4) = (A2*B2*(X(4)-Y02)**(B2-1.0))
  398. JACOB(6,4) = A2*B2*(A2*(X(4)-Y02))**(B2-1.0)
  399. JACOB(6,5) = 0.D0
  400. * DELAFON(6) = X(6)-(A2*(X(4)-Y02)**B2+1.0)
  401. DELAFON(6) = X(6)-((A2*(X(4)-Y02))**B2+1.0)
  402. DIMPB=6
  403. ELSE
  404. DIMPB=5
  405. END IF
  406. ELSE
  407. C
  408. C =====> PAS D'ENDOMMAGEMENT EN TRACTION
  409. C
  410. IF (LDOM2) THEN
  411. C
  412. C =====> ENDOMMAGEMENT EN COMPRESSION
  413. C
  414. J16 = (BETA2*NU+BETA1*(H2+X(5)*(1-H2)-1)*HT2+BETA1*(H2+X(5)*(1-H2)
  415. 1 -1)*HT+BETA2*(X(5)*H2+H2+2*X(6)*(1-H2)-1)+E*EPSP1*(H2-1))*SIGF1/
  416. 2 ((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X(5)**2*H2+X(6)**2*(1-H2))-X(5)*
  417. 3 X(6)*H2+X(6)**2*(H2-1))*SIGF1+(2-2*X(5))*BETA1*HT2*NU+(BETA1*(-X
  418. 4 (5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-2)+2))+((X(5)-1)*X(6)-X(5)
  419. 5 **2+X(5))*BETA1*H1)*HT2)-((H1*(X(5)*(2*H2-1)+2*X(6)*(1-H2))-X(5
  420. 6 )*H2+2*X(6)*(H2-1))*SIGF1+(BETA1*(-H2+X(5)*(H2-2)+2)+(X(5)-1)*BE
  421. 7 TA1*H1)*HT2)*((((X(5)-1)*BETA1*HT2+(X(5)-1)*BETA1*HT-E*EPSP2+(X
  422. 8 (6)-1)*BETA2)*NU+BETA1*(X(6)*(H2+X(5)*(1-H2)-1)+X(5)**2*H2-X(5)
  423. 9 *H2)*HT2+BETA1*(X(6)*(H2+X(5)*(1-H2)-1)+X(5)**2*H2-X(5)*H2)*HT+
  424. : BETA2*(X(6)*(X(5)*H2+H2-1)-X(5)*H2+X(6)**2*(1-H2))+E*EPSP1*(X(6
  425. ; )*(H2-1)-X(5)*H2))*SIGF1+(BETA1*E*(X(5)*EPSP2-EPSP2)+(1-X(5))*BE
  426. < TA1*E*EPSP1)*HT2)/((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X(5)**2*H2+X(6
  427. = )**2*(1-H2))-X(5)*X(6)*H2+X(6)**2*(H2-1))*SIGF1+(2-2*X(5))*BETA1
  428. > *HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-2)+2))+(
  429. ? (X(5)-1)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)**2
  430. C
  431. J26 = (BETA2*NU+((1-X(5))*BETA1*H1+(X(5)-1)*BETA1)*HT2+((1-X(5))*B
  432. 1 ETA1*H1+(X(5)-1)*BETA1)*HT+(E*EPSP2+(-2*X(6)+X(5)+1)*BETA2)*H1-
  433. 2 E*EPSP2+(2*X(6)-1)*BETA2)*SIGF1/((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X
  434. 3 (5)**2*H2+X(6)**2*(1-H2))-X(5)*X(6)*H2+X(6)**2*(H2-1))*SIGF1+(2-
  435. 4 2*X(5))*BETA1*HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)
  436. 5 *(H2-2)+2))+((X(5)-1)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)-((H1*(X
  437. 6 (5)*(2*H2-1)+2*X(6)*(1-H2))-X(5)*H2+2*X(6)*(H2-1))*SIGF1+(BETA1*
  438. 7 (-H2+X(5)*(H2-2)+2)+(X(5)-1)*BETA1*H1)*HT2)*((((X(5)-1)*BETA1*H
  439. 8 T2+(X(5)-1)*BETA1*HT-E*EPSP1+(X(6)-1)*BETA2)*NU+(((1-X(5))*X(6)
  440. 9 +X(5)**2-X(5))*BETA1*H1+(X(5)-1)*X(6)*BETA1)*HT2+(((1-X(5))*X(6
  441. : )+X(5)**2-X(5))*BETA1*H1+(X(5)-1)*X(6)*BETA1)*HT+(E*(X(6)*EPSP2
  442. ; -X(5)*EPSP2)+(-X(6)**2+(X(5)+1)*X(6)-X(5))*BETA2)*H1-X(6)*E*EPS
  443. < P2+(X(6)**2-X(6))*BETA2)*SIGF1+(BETA1*E*(EPSP2-X(5)*EPSP2)+(X(5)
  444. = -1)*BETA1*E*EPSP1)*HT2)/((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X(5)**2*
  445. > H2+X(6)**2*(1-H2))-X(5)*X(6)*H2+X(6)**2*(H2-1))*SIGF1+(2-2*X(5))
  446. ? *BETA1*HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-2)
  447. @ +2))+((X(5)-1)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)**2
  448. C
  449. C
  450.  
  451. JACOB(1,5) = J16
  452. JACOB(2,5) = J26
  453. JACOB(3,5) = 0.D0
  454. JACOB(4,5) = ((1.D0-H1)*X(1)*X(1)+(1.D0-H2)
  455. & *X(2)*X(2))*X(6)/E
  456. JACOB(5,5) = -1.D0
  457. & -A2*B2*(A2*(X(4)-Y02))**(B2-UN)
  458. & *alpha2*m2*dy02**alpha2/(dy02+hdp2)/(x(6)**2*dt1)
  459. JACOB(5,1) = 0.D0
  460. JACOB(5,2) = 0.D0
  461. JACOB(5,3) = 0.D0
  462. * JACOB(5,4) = (A2*B2*(X(4)-Y02)**(B2-1.0))
  463. JACOB(5,4) = A2*B2*(A2*(X(4)-Y02))**(B2-1.0)
  464. * DELAFON(5) = X(6)-(A2*(X(4)-Y02)**B2+1.0)
  465. DELAFON(5) = X(6)-((A2*(X(4)-Y02))**B2+1.0)
  466. DIMPB=5
  467. ELSE
  468. DIMPB=4
  469. END IF
  470. END IF
  471. C
  472. C RESOLUTION
  473. C
  474. CALL GAUSCL(JACOB,DX,DELAFON,NDIM,DIMPB,0)
  475. C
  476. C CALCUL DE L"ERREUR ET DE LA NOUVELLE SOLUTION
  477. C
  478. ERREUR=0.D0
  479. do i=1,ndim
  480. erreur=erreur+(delafon(i)*delafon(i))
  481. enddo
  482. ERREUR=0.D0
  483. DO I=1,2
  484. ERREUR=ERREUR+DX(I)*DX(I)/(SIGY*SIGY)
  485. X(I)=X(I)+DX(I)
  486. c if (x(i) .ge. 0.) then
  487. c x(i)=MIN(1.d+15,x(i))
  488. c else
  489. c x(i)= MAX(-1.d+15,x(i))
  490. c endif
  491. END DO
  492. X(3)=X(3)+DX(3)
  493. X3= MAX(X(3),Y01)
  494. c X3=MIN(ABS(X3),1.d+07)
  495. ERREUR=ERREUR+DX(3)*DX(3)/(x3*x3)
  496. X(4)=X(4)+DX(4)
  497. X4= MAX(X(4),Y02)
  498. c X4=MIN(ABS(X4),1.d+07)
  499. ERREUR=ERREUR+DX(4)*DX(4)/(x4*x4)
  500. c if (.not. ldom1) then
  501. c dx5=(x(5)-x5)*1.0/4.0
  502. c x(5)=dx5+x5
  503. c x(5)= MAX(z1,x(5))
  504. c ERREUR=ERREUR+DX5*DX5/(z1*z1)
  505. c endif
  506. c if (.not. ldom2) then
  507. c dx6=(x(6)-x6)/4.
  508. c x(6)=dx6+x6
  509. c x(6)= MAX(z2,x(6))
  510. c ERREUR=ERREUR+DX6*DX6/(z2*z2)
  511. c endif
  512. IF (LDOM1) THEN
  513. IF ((X(5)+DX(5)).LT.Z1) THEN
  514. DX(5)=(Z1-X(5))
  515. END IF
  516. IF ( X(5) .GT. ZMAX)THEN
  517. DX(5)=ZMAX + 0.1D0 - X(5)
  518. END IF
  519. X(5)=X(5)+DX(5)
  520. ERREUR=ERREUR+DX(5)*DX(5)/(z1*z1)
  521. IF (LDOM2) THEN
  522. IF ((X(6)+DX(6)).LT.Z2) THEN
  523. DX(6)=(Z2-X(6))
  524. END IF
  525. IF ( X(6) .GT. ZMAX)THEN
  526. DX(6)=ZMAX + 0.1D0 - X(6)
  527. END IF
  528. X(6)=X(6)+DX(6)
  529. ERREUR=ERREUR+DX(6)*DX(6)/(z2*z2)
  530. END IF
  531. ELSE IF (LDOM2) THEN
  532. IF ((X(6)+DX(5)).LT.Z2) THEN
  533. DX(5)=(Z2-X(6))
  534. END IF
  535. IF ( X(6) .GT. ZMAX)THEN
  536. DX(5)=ZMAX + 0.1D0 - X(6)
  537. END IF
  538. X(6)=X(6)+DX(5)
  539. ERREUR=ERREUR+DX(5)*DX(5)/(z2*z2)
  540. END IF
  541.  
  542.  
  543. IF (ERREUR .GT. 1.D+100)THEN
  544. CONV=.FALSE.
  545. print*,'erreur dans modele',erreur
  546. print*,'iteration num',iter
  547. print*,'delafon',delafon
  548. RETURN
  549. END IF
  550. C
  551. C QUELQUES IMPRESSIONS EN CAS D'ERREUR
  552. C
  553. if (.not. conv )then
  554. PRINT*,'ITERATION ',ITER,'ERREUR',erreur
  555. print*,'LDOM1 :',LDOM1,' LDOM2 :',LDOM2
  556. PRINT*,'EPS1= ',EPSP1,' EPS2= ',EPSP2
  557. PRINT*,'SIG1:',X(1) ,' DSIG1:',DX(1),'ERR :',DELAFON(1)
  558. PRINT*,'SIG2:',X(2) ,' DSIG2:',DX(2),'ERR :',DELAFON(2)
  559. PRINT*,'Y1 :',X(3) ,' DY1 :',DX(3),'ERR :',DELAFON(3)
  560. PRINT*,'Y2 :',X(4) ,' DY2 :',DX(4),'ERR :',DELAFON(4)
  561. PRINT*,'Z1 :',X(5) ,' DZ1 :',DX(5),'ERR :',DELAFON(5)
  562. PRINT*,'Z2 :',X(6) ,' DZ2 :',DX(6),'ERR :',DELAFON(6)
  563. PRINT*,'H1',H1,'H2',H2,ht,ht2
  564. PRINT*
  565. write(*,200)jacob
  566. END IF
  567.  
  568. END DO
  569.  
  570. *** pour z2 ****
  571. ********************
  572. fctd2=0.
  573. If ( x(4) .Ge. Y02 ) Then
  574. ** si z converge est > a x(6) alors fctd est < 0
  575. Fctd2 =-( 1.0+(A2*(-Y02 + X(4)))**B2-X(6) )
  576. if (fctd2.ge.0.0 .and. x(6).eq.z2) fctd2 =0.0
  577. Else
  578. if (x(4) .le. y002 ) then
  579. fctd2=x(6)-z2
  580. c x(6)=z2
  581. else
  582. Fctd2 =x(6)-z2
  583. endif
  584. Endif
  585. If ( Fctd2 .lt. 0.0) then
  586. Xg2=Xd2
  587. if (passage2) pas2=pas2*2.0
  588. Xd2=Xd2+pas2
  589. Else
  590. Xd2=Xg2
  591. Pas2=Pas2/2.0
  592. Xd2=Xd2+pas2
  593. passage2=.false.
  594. Endif
  595. Erreurd2=(fctd2*fctd2)/(z2*z2)
  596. if (x(6) .gt. zmax) then
  597. x(6)=zmax
  598. Erreurd2=0.0
  599. endif
  600.  
  601. c print*,'Erreurd2'
  602. c print*,Erreurd2,x(4),y02,y002
  603.  
  604. ENDDO
  605.  
  606. *** pour z1 ****
  607. ********************
  608. fctd1=0.0
  609. If ( x(3) .Ge. Y01 ) Then
  610. ** si z converge est > a x(5) alors fctd1 est < 0
  611. Fctd1 =-( 1.0+(A1*(-Y01 + X(3)))**B1-X(5) )
  612. if (fctd1.ge.0.0 .and. x(5).eq.z1) fctd1=0.0
  613. Else
  614. if (x(3) .le. y001) then
  615. fctd1=x(5)-z2
  616. x(5)=z1
  617. else
  618. Fctd1 =x(5)-z2
  619. endif
  620. Endif
  621. If ( Fctd1 .lt. 0.0) then
  622. Xg1=Xd1
  623. if (passage1) pas1=pas1*2.0
  624. Xd1=Xd1+pas1
  625. Else
  626. Xd1=Xg1
  627. Pas1=Pas1/2.0
  628. Xd1=Xd1+pas1
  629. passage1=.false.
  630. Endif
  631. Erreurd1=(Fctd1*Fctd1)/(z1*z1)
  632. if (x(5) .gt. zmax) then
  633. x(5)=zmax
  634. Erreurd1=0.0
  635. endif
  636.  
  637. c print*,'Erreurd1'
  638. c print*,Erreurd1,x(3),y02,Y002
  639.  
  640. ENDDO
  641.  
  642. IF (ERREUR .GT. TEST ) THEN
  643. CONV=.FALSE.
  644. ELSE
  645. CONV=.TRUE.
  646. SIG1= X(1)
  647. SIG2= X(2)
  648. Y1 = X(3)
  649. Y2 = X(4)
  650. Z1 = X(5)
  651. Z2 = X(6)
  652. IF (SIG1 .GE. 0.D0) THEN
  653. JSIGN(1)=1
  654. ELSE
  655. JSIGN(1)=0
  656. END IF
  657. IF (SIG2.GE. 0.D0)THEN
  658. JSIGN(2)=1
  659. ELSE
  660. JSIGN(2)=0
  661. END IF
  662. JSIGN(3)=1
  663. END IF
  664. RETURN
  665. 200 format(6(6(2x,e11.4),/))
  666. END
  667.  
  668.  
  669.  
  670.  
  671.  
  672.  
  673.  
  674.  
  675.  
  676.  
  677.  
  678.  
  679.  
  680.  
  681.  
  682.  

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