Télécharger codjfd.eso

Retour à la liste

Numérotation des lignes :

  1. C CODJFD SOURCE BP208322 17/03/01 21:15:59 9325
  2. C MODJFD SOURCE AM 00/12/13 21:41:28 4045
  3. SUBROUTINE CODJFD(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 MODJFD 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 NEWTON %
  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. UN=1.D0
  61. DEUX=2.D0
  62. ZMAX=1.d+04
  63. E = XMAT(1)
  64. NU = XMAT(2)
  65. Y001 = XMAT(5)
  66. Y002 = XMAT(6)
  67. A1 = XMAT(7)
  68. A2 = XMAT(8)
  69. B1 = XMAT(9)
  70. B2 = XMAT(10)
  71. BETA1 = XMAT(11)
  72. BETA2 = XMAT(12)
  73. SIGF1 = XMAT(13)
  74. M1 = XMAT(14)
  75. M2 = XMAT(15)
  76. ALPHA1= XMAT(16)
  77. ALPHA2= XMAT(17)
  78. SIGY =SQRT(2.D0*E*Y001)
  79. * DELTASIG=SIGY*1.D-4
  80. DELTASIG=0.d0
  81. * print*,'E,NU,Y01,Y02,A1,A2,B1,B2,BETA1,BETA2'
  82. * print*,E,NU,Y01,Y02,A1,A2,B1,B2,BETA1,BETA2
  83. TEST=1.D-12
  84. ERREUR=1.D0
  85. ITER=0
  86. X(1)=SIG1
  87. X(2)=SIG2
  88. X(3)=Y1
  89. X(4)=Y2
  90. X(5)=Z1
  91. X(6)=Z2
  92. DO I=1, NDIM
  93. X0(I)=X(I)
  94. END DO
  95. x5=Z1
  96. x6=z2
  97. C
  98. C
  99. C
  100. JACOB(1,1) = -1.d0
  101. JACOB(1,2) = 0.D0
  102. JACOB(1,3) = 0.D0
  103. JACOB(1,4) = 0.D0
  104. JACOB(1,5) = 0.D0
  105. JACOB(1,6) = 0.D0
  106. JACOB(2,1) = 0.D0
  107. JACOB(2,2) = -1.d0
  108. JACOB(2,3) = 0.D0
  109. JACOB(2,4) = 0.D0
  110. JACOB(2,5) = 0.D0
  111. JACOB(2,6) = 0.D0
  112. JACOB(3,1) = 0.D0
  113. JACOB(3,2) = 0.D0
  114. JACOB(3,3) = -1.d0
  115. JACOB(3,4) = 0.D0
  116. JACOB(3,5) = 0.D0
  117. JACOB(3,6) = 0.D0
  118. JACOB(4,1) = 0.D0
  119. JACOB(4,2) = 0.D0
  120. JACOB(4,3) = 0.D0
  121. JACOB(4,4) = -1.d0
  122. JACOB(4,5) = 0.D0
  123. JACOB(4,6) = 0.D0
  124. JACOB(5,1) = 0.D0
  125. JACOB(5,2) = 0.D0
  126. JACOB(5,3) = 0.D0
  127. JACOB(5,4) = 0.D0
  128. JACOB(5,5) = -1.d0
  129. JACOB(5,6) = 0.D0
  130. JACOB(6,1) = 0.D0
  131. JACOB(6,2) = 0.D0
  132. JACOB(6,3) = 0.D0
  133. JACOB(6,4) = 0.D0
  134. JACOB(6,5) = 0.D0
  135. JACOB(6,6) = -1.d0
  136.  
  137. DO WHILE (ERREUR .GT. TEST .AND. ITER .LE. 100)
  138. ITER=ITER+1
  139. IF (X(1) .GE. DELTASIG) THEN
  140. H1=UN
  141. ELSE IF (X(1) .LE. -DELTASIG) THEN
  142. H1=0.D0
  143. ELSE
  144. H1=(X(1)+DELTASIG)/(DEUX*DELTASIG)
  145. END IF
  146. IF (X(2) .GE. DELTASIG) THEN
  147. H2=UN
  148. ELSE IF (X(2) .LE. -DELTASIG) THEN
  149. H2=0.D0
  150. ELSE
  151. H2=(X(2)+DELTASIG)/(DEUX*DELTASIG)
  152. END IF
  153. c
  154. c smooth transition on ht ht2
  155. c
  156. trf = x(1) + x(2)
  157. IF(trf.GE.deltasig)THEN
  158. HT=UN
  159. HT2=0.D0
  160. else if(trf.ge.(-deltasig)) then
  161. ht = (trf+deltasig)/DEUX/deltasig
  162. ht2 = UN-ht
  163. else if(trf.ge.(-sigf+deltasig)) then
  164. ht = 0.D+00
  165. ht2 = UN
  166. else if(trf.ge.(-sigf-deltasig)) then
  167. ht = 0.D+00
  168. c ht2 =UN- (trf-(sigf-deltasig))/(DEUX*sigf)
  169. ht2 = (trf+sigf+deltasig)/DEUX/deltasig
  170. else
  171. ht = 0.D+00
  172. ht2 = 0.D+00
  173. endif
  174. c print*,'DTjfd',dt1
  175. dy01=(ddim(x(5),z1)/(x(5)*z1*dt1))
  176. hdp1=1.0
  177. if (dy01 .gt. 0.0) hdp1=0.0
  178. y01=y001+m1*dy01**alpha1
  179. dy02=(ddim(x(6),z2)/(x(6)*z2*dt1))
  180. hdp2=1.0
  181. if (dy02 .gt. 0.0) hdp2=0.0
  182. y02=y002+m2*dy02**alpha2
  183. c print*,'iter=',iter,' dy01=',dy01,'dy02=',dy02
  184.  
  185. ldom1=.false.
  186. evol1=.false.
  187. if (X(3).GT. Y01 .AND. X(5) .LT. ZMAX) then
  188. if ( (1.+(A1*(-Y01 + X(3)))**B1-z1) .gt. 1.d-06) then
  189. evol1=.true.
  190. ldom1=.true.
  191. x5=z1
  192. endif
  193. endif
  194. ldom2=.false.
  195. evol2=.false.
  196. if (X(4).GT. Y02 .AND. X(6) .LT. ZMAX) then
  197. if ( (1.+(A2*(-Y02 + X(4)))**B2-z2) .gt. 1.d-06) then
  198. evol2=.true.
  199. ldom2=.true.
  200. x6=z2
  201. endif
  202. endif
  203. c print*,'y1',x(3),'ylim1',ylim1
  204. c print*,'z1',x(5),'x5',x5
  205. c print*,'ldom1',ldom1,'ldom2',ldom2
  206. C =====> PROBLEME ELASTIQUE
  207. C
  208. C
  209. C CALCUL DU JACOBIEN
  210. C
  211. C
  212. JACOB(3,1) = H1*X(1)*X(5)*X(5)/E
  213. C
  214. JACOB(3,2) = H2*X(2)*X(5)*X(5)/E
  215. C
  216. JACOB(4,1) = (UN-H1)*X(1)*X(6)*X(6)/E
  217. C
  218. JACOB(4,2) = (UN-H2)*X(2)*X(6)*X(6)/E
  219. C
  220. C CALCUL DU SECOND MEMBRE
  221. C
  222. XIG1 =((((X(5)-UN)*BETA1*HT2+(X(5)-UN)*BETA1*HT-E*EPSP2+(X(6)-UN)
  223. 1 *BETA2)*NU+BETA1*(X(6)*(H2+X(5)*(UN-H2)-UN)+X(5)**2*H2-X(5)*H2)
  224. 2 *HT2+BETA1*(X(6)*(H2+X(5)*(UN-H2)-UN)+X(5)**2*H2-X(5)*H2)*HT+
  225. 3 BETA2*(X(6)*(X(5)*H2+H2-UN)-X(5)*H2+X(6)**2*(UN-H2))+E*EPSP1*
  226. 4 (X(6)*(H2-UN)-X(5)*H2))*SIGF1+(BETA1*E*(X(5)*EPSP2-EPSP2)+
  227. 5 (UN-X(5))*BETA1*E*EPSP1)*HT2)/((NU**2+H1*(X(5)*X(6)*(2*H2-UN)-
  228. 6 X(5)**2*H2+X(6)**2*(UN-H2))-X(5)*X(6)*H2+X(6)**2*(H2-UN))*SIGF1+
  229. 7 (DEUX-DEUX*X(5))*BETA1*HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*
  230. 8 (-H2+X(5)*(H2-DEUX)+DEUX))+((X(5)-UN)*X(6)-X(5)**2+X(5))
  231. 9 *BETA1*H1)*HT2)
  232. C
  233. XIG2 =((((X(5)-UN)*BETA1*HT2+(X(5)-UN)*BETA1*HT-E*EPSP1+(X(6)-UN)
  234. 1 *BETA2)*NU+(((UN-X(5))*X(6)+X(5)**2-X(5))*BETA1*H1+(X(5)-UN)*
  235. 2 X(6)*BETA1)*HT2+(((UN-X(5))*X(6)+X(5)**2-X(5))*BETA1*H1+
  236. 3 (X(5)-UN)*X(6)*BETA1)*HT+(E*(X(6)*EPSP2-X(5)*EPSP2)+(-X(6)**2+
  237. 4 (X(5)+UN)*X(6)-X(5))*BETA2)*H1-X(6)*E*EPSP2+(X(6)**2-X(6))*
  238. 5BETA2)*SIGF1+(BETA1*E*(EPSP2-X(5)*EPSP2)+(X(5)-UN)*BETA1*E*EPSP1)
  239. 6 *HT2)/((NU**2+H1*(X(5)*X(6)*(DEUX*H2-UN)-X(5)**2*H2+X(6)**2*
  240. 7 (UN-H2))-X(5)*X(6)*H2+X(6)**2*(H2-UN))*SIGF1+(DEUX-DEUX*X(5))*
  241. 8 BETA1*HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*
  242. 9 (H2-DEUX)+DEUX))+((X(5)-UN)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)
  243.  
  244. DELAFON(1) =X(1)-XIG1
  245. C
  246. DELAFON(2) =X(2)-XIG2
  247. C
  248. XX= (H1*X(1)*X(1)+H2*X(2)*X(2))*X(5)*X(5)/(DEUX*E)
  249. DELAFON(3) = X(3)-XX
  250. C
  251. XX=((UN-H1)*X(1)*X(1)+(UN-H2)*X(2)*X(2) )
  252. 1 *X(6)*X(6)/(DEUX*E)
  253. DELAFON(4) = X(4)-XX
  254. C
  255. IF (LDOM1) THEN
  256. C
  257. C =====> ENDOMMAGEMENT EN TRACTION
  258. C
  259. J15 = (((BETA1*HT2+BETA1*HT)*NU+BETA1*(DEUX*X(5)*H2-H2+X(6)*(UN-H2
  260. & ))*HT2+BETA1*(DEUX*X(5)*H2-H2+X(6)*(UN-H2))*HT+BETA2*(X(6)*H2-H
  261. & 2)-E*EPSP1*H2)*SIGF1+(BETA1*E*EPSP2-BETA1*E*EPSP1)*HT2)/((NU**2+
  262. & H1*(X(5)*X(6)*(DEUX*H2-UN)-X(5)**2*H2+X(6)**2*(UN-H2))-X(5)*X(6
  263. & )*H2+X(6)**2*(H2-UN))*SIGF1+(DEUX-DEUX*X(5))*BETA1*HT2*NU+(BETA1
  264. & *(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-DEUX)+DEUX))+((X(5)-UN
  265. & )*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)-((H1*(X(6)*(DEUX*H2-UN)-DEU
  266. & X*X(5)*H2)-X(6)*H2)*SIGF1-DEUX*BETA1*HT2*NU+(BETA1*(-DEUX*X(5)*H
  267. & 2+H2+X(6)*(H2-DEUX))+(X(6)-DEUX*X(5)+UN)*BETA1*H1)*HT2)*((((X(5
  268. & )-UN)*BETA1*HT2+(X(5)-UN)*BETA1*HT-E*EPSP2+(X(6)-UN)*BETA2)*NU+
  269. & BETA1*(X(6)*(H2+X(5)*(UN-H2)-UN)+X(5)**2*H2-X(5)*H2)*HT2+BETA1*
  270. & (X(6)*(H2+X(5)*(UN-H2)-UN)+X(5)**2*H2-X(5)*H2)*HT+BETA2*(X(6)*(
  271. & X(5)*H2+H2-UN)-X(5)*H2+X(6)**2*(UN-H2))+E*EPSP1*(X(6)*(H2-UN)-X
  272. & (5)*H2))*SIGF1+(BETA1*E*(X(5)*EPSP2-EPSP2)+(UN-X(5))*BETA1*E*EPS
  273. & P1)*HT2)/((NU**2+H1*(X(5)*X(6)*(DEUX*H2-UN)-X(5)**2*H2+X(6)**2*
  274. & (UN-H2))-X(5)*X(6)*H2+X(6)**2*(H2-UN))*SIGF1+(DEUX-DEUX*X(5))*BE
  275. & TA1*HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-DEUX)
  276. & +DEUX))+((X(5)-UN)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)**2
  277. C
  278. C
  279. J25 = (((BETA1*HT2+BETA1*HT)*NU+((-X(6)+DEUX*X(5)-UN)*BETA1*H1+X(6
  280. & )*BETA1)*HT2+((-X(6)+DEUX*X(5)-UN)*BETA1*H1+X(6)*BETA1)*HT+((X(
  281. & 6)-UN)*BETA2-E*EPSP2)*H1)*SIGF1+(BETA1*E*EPSP1-BETA1*E*EPSP2)*HT
  282. & 2)/((NU**2+H1*(X(5)*X(6)*(DEUX*H2-UN)-X(5)**2*H2+X(6)**2*(UN-H2
  283. & ))-X(5)*X(6)*H2+X(6)**2*(H2-UN))*SIGF1+(DEUX-DEUX*X(5))*BETA1*HT
  284. & 2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-DEUX)+DEUX)
  285. & )+((X(5)-UN)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)-((H1*(X(6)*(DEUX
  286. & *H2-UN)-DEUX*X(5)*H2)-X(6)*H2)*SIGF1-DEUX*BETA1*HT2*NU+(BETA1*(-
  287. & DEUX*X(5)*H2+H2+X(6)*(H2-DEUX))+(X(6)-DEUX*X(5)+UN)*BETA1*H1)*H
  288. & T2)*((((X(5)-UN)*BETA1*HT2+(X(5)-UN)*BETA1*HT-E*EPSP1+(X(6)-UN)
  289. & *BETA2)*NU+(((UN-X(5))*X(6)+X(5)**2-X(5))*BETA1*H1+(X(5)-UN)*X(
  290. & 6)*BETA1)*HT2+(((UN-X(5))*X(6)+X(5)**2-X(5))*BETA1*H1+(X(5)-UN)
  291. & *X(6)*BETA1)*HT+(E*(X(6)*EPSP2-X(5)*EPSP2)+(-X(6)**2+(X(5)+UN)*
  292. & X(6)-X(5))*BETA2)*H1-X(6)*E*EPSP2+(X(6)**2-X(6))*BETA2)*SIGF1+(B
  293. & ETA1*E*(EPSP2-X(5)*EPSP2)+(X(5)-UN)*BETA1*E*EPSP1)*HT2)/((NU**2
  294. & +H1*(X(5)*X(6)*(DEUX*H2-UN)-X(5)**2*H2+X(6)**2*(UN-H2))-X(5)*X(
  295. & 6)*H2+X(6)**2*(H2-UN))*SIGF1+(DEUX-DEUX*X(5))*BETA1*HT2*NU+(BETA
  296. & 1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-DEUX)+DEUX))+((X(5)-U
  297. & N)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)**2
  298. C
  299. JACOB(1,5) = J15
  300. JACOB(2,5) = J25
  301. JACOB(3,5) = (H1*X(1)*X(1)+H2*X(2)*X(2))*X(5)/E
  302. JACOB(4,5) = 0.D0
  303. JACOB(5,5) = -1.D0
  304. & -A1*B1*(A1*(X(3)-Y01))**(B1-UN)
  305. & *alpha1*m1*dy01**alpha1/(dy01+hdp1)/(x(5)**2*dt1)
  306. JACOB(5,1) = 0.D0
  307. JACOB(5,2) = 0.D0
  308. JACOB(5,3) = A1*B1*(A1*(X(3)-Y01))**(B1-UN)
  309. JACOB(5,4) = 0.D0
  310. DELAFON(5) = X(5)-((A1*(X(3)-Y01))**B1+UN)
  311.  
  312. IF (LDOM2) THEN
  313. C
  314. C =====> ENDOMMAGEMENT EN COMPRESSION
  315. C
  316. J16 = (BETA2*NU+BETA1*(H2+X(5)*(1-H2)-1)*HT2+BETA1*(H2+X(5)*(1-H2)
  317. 1 -1)*HT+BETA2*(X(5)*H2+H2+2*X(6)*(1-H2)-1)+E*EPSP1*(H2-1))*SIGF1/
  318. 2 ((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X(5)**2*H2+X(6)**2*(1-H2))-X(5)*
  319. 3 X(6)*H2+X(6)**2*(H2-1))*SIGF1+(2-2*X(5))*BETA1*HT2*NU+(BETA1*(-X
  320. 4 (5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-2)+2))+((X(5)-1)*X(6)-X(5)
  321. 5 **2+X(5))*BETA1*H1)*HT2)-((H1*(X(5)*(2*H2-1)+2*X(6)*(1-H2))-X(5
  322. 6 )*H2+2*X(6)*(H2-1))*SIGF1+(BETA1*(-H2+X(5)*(H2-2)+2)+(X(5)-1)*BE
  323. 7 TA1*H1)*HT2)*((((X(5)-1)*BETA1*HT2+(X(5)-1)*BETA1*HT-E*EPSP2+(X
  324. 8 (6)-1)*BETA2)*NU+BETA1*(X(6)*(H2+X(5)*(1-H2)-1)+X(5)**2*H2-X(5)
  325. 9 *H2)*HT2+BETA1*(X(6)*(H2+X(5)*(1-H2)-1)+X(5)**2*H2-X(5)*H2)*HT+
  326. : BETA2*(X(6)*(X(5)*H2+H2-1)-X(5)*H2+X(6)**2*(1-H2))+E*EPSP1*(X(6
  327. ; )*(H2-1)-X(5)*H2))*SIGF1+(BETA1*E*(X(5)*EPSP2-EPSP2)+(1-X(5))*BE
  328. < TA1*E*EPSP1)*HT2)/((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X(5)**2*H2+X(6
  329. = )**2*(1-H2))-X(5)*X(6)*H2+X(6)**2*(H2-1))*SIGF1+(2-2*X(5))*BETA1
  330. > *HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-2)+2))+(
  331. ? (X(5)-1)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)**2
  332. C
  333. J26 = (BETA2*NU+((1-X(5))*BETA1*H1+(X(5)-1)*BETA1)*HT2+((1-X(5))*B
  334. 1 ETA1*H1+(X(5)-1)*BETA1)*HT+(E*EPSP2+(-2*X(6)+X(5)+1)*BETA2)*H1-
  335. 2 E*EPSP2+(2*X(6)-1)*BETA2)*SIGF1/((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X
  336. 3 (5)**2*H2+X(6)**2*(1-H2))-X(5)*X(6)*H2+X(6)**2*(H2-1))*SIGF1+(2-
  337. 4 2*X(5))*BETA1*HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)
  338. 5 *(H2-2)+2))+((X(5)-1)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)-((H1*(X
  339. 6 (5)*(2*H2-1)+2*X(6)*(1-H2))-X(5)*H2+2*X(6)*(H2-1))*SIGF1+(BETA1*
  340. 7 (-H2+X(5)*(H2-2)+2)+(X(5)-1)*BETA1*H1)*HT2)*((((X(5)-1)*BETA1*H
  341. 8 T2+(X(5)-1)*BETA1*HT-E*EPSP1+(X(6)-1)*BETA2)*NU+(((1-X(5))*X(6)
  342. 9 +X(5)**2-X(5))*BETA1*H1+(X(5)-1)*X(6)*BETA1)*HT2+(((1-X(5))*X(6
  343. : )+X(5)**2-X(5))*BETA1*H1+(X(5)-1)*X(6)*BETA1)*HT+(E*(X(6)*EPSP2
  344. ; -X(5)*EPSP2)+(-X(6)**2+(X(5)+1)*X(6)-X(5))*BETA2)*H1-X(6)*E*EPS
  345. < P2+(X(6)**2-X(6))*BETA2)*SIGF1+(BETA1*E*(EPSP2-X(5)*EPSP2)+(X(5)
  346. = -1)*BETA1*E*EPSP1)*HT2)/((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X(5)**2*
  347. > H2+X(6)**2*(1-H2))-X(5)*X(6)*H2+X(6)**2*(H2-1))*SIGF1+(2-2*X(5))
  348. ? *BETA1*HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-2)
  349. @ +2))+((X(5)-1)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)**2
  350. C
  351. C
  352. JACOB(1,6) = J16
  353. JACOB(2,6) = J26
  354. JACOB(3,6) = 0.D0
  355. JACOB(4,6) = ((1.D0-H1)*X(1)*X(1)+(1.D0-H2)
  356. & *X(2)*X(2))*X(6)/E
  357. JACOB(5,6) = 0.D0
  358. JACOB(6,6) = -1.D0
  359. & -A2*B2*(A2*(X(4)-Y02))**(B2-UN)
  360. & *alpha2*m2*dy02**alpha2/(dy02+hdp2)/(x(6)**2*dt1)
  361. JACOB(6,1) = 0.D0
  362. JACOB(6,2) = 0.D0
  363. JACOB(6,3) = 0.D0
  364. * JACOB(6,4) = (A2*B2*(X(4)-Y02)**(B2-1.0))
  365. JACOB(6,4) = A2*B2*(A2*(X(4)-Y02))**(B2-1.0)
  366. JACOB(6,5) = 0.D0
  367. * DELAFON(6) = X(6)-(A2*(X(4)-Y02)**B2+1.0)
  368. DELAFON(6) = X(6)-((A2*(X(4)-Y02))**B2+1.0)
  369. DIMPB=6
  370. ELSE
  371. DIMPB=5
  372. END IF
  373. ELSE
  374. C
  375. C =====> PAS D'ENDOMMAGEMENT EN TRACTION
  376. C
  377. IF (LDOM2) THEN
  378. C
  379. C =====> ENDOMMAGEMENT EN COMPRESSION
  380. C
  381. J16 = (BETA2*NU+BETA1*(H2+X(5)*(1-H2)-1)*HT2+BETA1*(H2+X(5)*(1-H2)
  382. 1 -1)*HT+BETA2*(X(5)*H2+H2+2*X(6)*(1-H2)-1)+E*EPSP1*(H2-1))*SIGF1/
  383. 2 ((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X(5)**2*H2+X(6)**2*(1-H2))-X(5)*
  384. 3 X(6)*H2+X(6)**2*(H2-1))*SIGF1+(2-2*X(5))*BETA1*HT2*NU+(BETA1*(-X
  385. 4 (5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-2)+2))+((X(5)-1)*X(6)-X(5)
  386. 5 **2+X(5))*BETA1*H1)*HT2)-((H1*(X(5)*(2*H2-1)+2*X(6)*(1-H2))-X(5
  387. 6 )*H2+2*X(6)*(H2-1))*SIGF1+(BETA1*(-H2+X(5)*(H2-2)+2)+(X(5)-1)*BE
  388. 7 TA1*H1)*HT2)*((((X(5)-1)*BETA1*HT2+(X(5)-1)*BETA1*HT-E*EPSP2+(X
  389. 8 (6)-1)*BETA2)*NU+BETA1*(X(6)*(H2+X(5)*(1-H2)-1)+X(5)**2*H2-X(5)
  390. 9 *H2)*HT2+BETA1*(X(6)*(H2+X(5)*(1-H2)-1)+X(5)**2*H2-X(5)*H2)*HT+
  391. : BETA2*(X(6)*(X(5)*H2+H2-1)-X(5)*H2+X(6)**2*(1-H2))+E*EPSP1*(X(6
  392. ; )*(H2-1)-X(5)*H2))*SIGF1+(BETA1*E*(X(5)*EPSP2-EPSP2)+(1-X(5))*BE
  393. < TA1*E*EPSP1)*HT2)/((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X(5)**2*H2+X(6
  394. = )**2*(1-H2))-X(5)*X(6)*H2+X(6)**2*(H2-1))*SIGF1+(2-2*X(5))*BETA1
  395. > *HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-2)+2))+(
  396. ? (X(5)-1)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)**2
  397. C
  398. J26 = (BETA2*NU+((1-X(5))*BETA1*H1+(X(5)-1)*BETA1)*HT2+((1-X(5))*B
  399. 1 ETA1*H1+(X(5)-1)*BETA1)*HT+(E*EPSP2+(-2*X(6)+X(5)+1)*BETA2)*H1-
  400. 2 E*EPSP2+(2*X(6)-1)*BETA2)*SIGF1/((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X
  401. 3 (5)**2*H2+X(6)**2*(1-H2))-X(5)*X(6)*H2+X(6)**2*(H2-1))*SIGF1+(2-
  402. 4 2*X(5))*BETA1*HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)
  403. 5 *(H2-2)+2))+((X(5)-1)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)-((H1*(X
  404. 6 (5)*(2*H2-1)+2*X(6)*(1-H2))-X(5)*H2+2*X(6)*(H2-1))*SIGF1+(BETA1*
  405. 7 (-H2+X(5)*(H2-2)+2)+(X(5)-1)*BETA1*H1)*HT2)*((((X(5)-1)*BETA1*H
  406. 8 T2+(X(5)-1)*BETA1*HT-E*EPSP1+(X(6)-1)*BETA2)*NU+(((1-X(5))*X(6)
  407. 9 +X(5)**2-X(5))*BETA1*H1+(X(5)-1)*X(6)*BETA1)*HT2+(((1-X(5))*X(6
  408. : )+X(5)**2-X(5))*BETA1*H1+(X(5)-1)*X(6)*BETA1)*HT+(E*(X(6)*EPSP2
  409. ; -X(5)*EPSP2)+(-X(6)**2+(X(5)+1)*X(6)-X(5))*BETA2)*H1-X(6)*E*EPS
  410. < P2+(X(6)**2-X(6))*BETA2)*SIGF1+(BETA1*E*(EPSP2-X(5)*EPSP2)+(X(5)
  411. = -1)*BETA1*E*EPSP1)*HT2)/((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X(5)**2*
  412. > H2+X(6)**2*(1-H2))-X(5)*X(6)*H2+X(6)**2*(H2-1))*SIGF1+(2-2*X(5))
  413. ? *BETA1*HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-2)
  414. @ +2))+((X(5)-1)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)**2
  415. C
  416. C
  417.  
  418. JACOB(1,5) = J16
  419. JACOB(2,5) = J26
  420. JACOB(3,5) = 0.D0
  421. JACOB(4,5) = ((1.D0-H1)*X(1)*X(1)+(1.D0-H2)
  422. & *X(2)*X(2))*X(6)/E
  423. JACOB(5,5) = -1.D0
  424. & -A2*B2*(A2*(X(4)-Y02))**(B2-UN)
  425. & *alpha2*m2*dy02**alpha2/(dy02+hdp2)/(x(6)**2*dt1)
  426. JACOB(5,1) = 0.D0
  427. JACOB(5,2) = 0.D0
  428. JACOB(5,3) = 0.D0
  429. * JACOB(5,4) = (A2*B2*(X(4)-Y02)**(B2-1.0))
  430. JACOB(5,4) = A2*B2*(A2*(X(4)-Y02))**(B2-1.0)
  431. * DELAFON(5) = X(6)-(A2*(X(4)-Y02)**B2+1.0)
  432. DELAFON(5) = X(6)-((A2*(X(4)-Y02))**B2+1.0)
  433. DIMPB=5
  434. ELSE
  435. DIMPB=4
  436. END IF
  437. END IF
  438. C
  439. C RESOLUTION
  440. C
  441. CALL GAUSCL(JACOB,DX,DELAFON,NDIM,DIMPB,0)
  442. C
  443. C CALCUL DE L"ERREUR ET DE LA NOUVELLE SOLUTION
  444. C
  445. ERREUR=0.D0
  446. DO I=1,2
  447. ERREUR=ERREUR+DX(I)*DX(I)/(SIGY*SIGY)
  448. X(I)=X(I)+DX(I)
  449. if (x(i) .ge. 0.) then
  450. x(i)=dmin1(1.d+15,x(i))
  451. else
  452. x(i)=dmax1(-1.d+15,x(i))
  453. endif
  454. END DO
  455. X(3)=X(3)+DX(3)
  456. X3=DMAX1(X(3),Y001)
  457. X3=Dmin1(ABS(X3),1.d+07)
  458. ERREUR=ERREUR+DX(3)*DX(3)/(x3*x3)
  459. X(4)=X(4)+DX(4)
  460. X4=DMAX1(X(4),Y002)
  461. X4=Dmin1(ABS(X4),1.d+07)
  462. ERREUR=ERREUR+DX(4)*DX(4)/(x4*x4)
  463. if (.not. ldom1) then
  464. dx5=(x(5)-x5)*1.0/4.0
  465. x(5)=dx5+x5
  466. x(5)=dmax1(z1,x(5))
  467. ERREUR=ERREUR+DX5*DX5/(z1*z1)
  468. endif
  469. if (.not. ldom2) then
  470. dx6=(x(6)-x6)/4.
  471. x(6)=dx6+x6
  472. x(6)=dmax1(z2,x(6))
  473. ERREUR=ERREUR+DX6*DX6/(z2*z2)
  474. endif
  475. IF (LDOM1) THEN
  476. IF ((X(5)+DX(5)).LT.Z1) THEN
  477. DX(5)=(Z1-X(5))
  478.  
  479. c if (evol1) dx(5)=dx(5)/2.
  480.  
  481. END IF
  482. IF ( X(5) .GT. ZMAX)THEN
  483. DX(5)=ZMAX + 0.1D0 - X(5)
  484. END IF
  485. X(5)=X(5)+DX(5)
  486. ERREUR=ERREUR+DX(5)*DX(5)/(z1*z1)
  487. IF (LDOM2) THEN
  488. IF ((X(6)+DX(6)).LT.Z2) THEN
  489. DX(6)=(Z2-X(6))
  490.  
  491. c if (evol2) dx(6)=dx(6)/2.
  492.  
  493. END IF
  494. IF ( X(6) .GT. ZMAX)THEN
  495. DX(6)=ZMAX + 0.1D0 - X(6)
  496. END IF
  497. X(6)=X(6)+DX(6)
  498. ERREUR=ERREUR+DX(6)*DX(6)/(z2*z2)
  499. END IF
  500. ELSE IF (LDOM2) THEN
  501. IF ((X(6)+DX(5)).LT.Z2) THEN
  502. DX(5)=(Z2-X(6))
  503.  
  504. c if (evol2) dx(5)=dx(5)/2.
  505.  
  506. END IF
  507. IF ( X(6) .GT. ZMAX)THEN
  508. DX(5)=ZMAX + 0.1D0 - X(6)
  509. END IF
  510. X(6)=X(6)+DX(5)
  511. ERREUR=ERREUR+DX(5)*DX(5)/(z2*z2)
  512. END IF
  513.  
  514.  
  515. IF (ERREUR .GT. 1.D+100)THEN
  516. CONV=.FALSE.
  517. print*,'erreur dans modele',erreur
  518. print*,'iteration num',iter
  519. print*,'delafon',delafon
  520. print*,z1,z2,x3,x4,sigy
  521. print*,dx
  522. print*,'y02=',y02,'y002=',y002
  523. RETURN
  524. END IF
  525. C
  526. C QUELQUES IMPRESSIONS EN CAS D'ERREUR
  527. C
  528. if (.not. conv )then
  529. PRINT*,'ITERATION ',ITER,'ERREUR',erreur
  530. print*,'LDOM1 :',LDOM1,' LDOM2 :',LDOM2
  531. PRINT*,'EPS1= ',EPSP1,' EPS2= ',EPSP2
  532. PRINT*,'SIG1:',X(1) ,' DSIG1:',DX(1),'ERR :',DELAFON(1)
  533. PRINT*,'SIG2:',X(2) ,' DSIG2:',DX(2),'ERR :',DELAFON(2)
  534. PRINT*,'Y1 :',X(3) ,' DY1 :',DX(3),'ERR :',DELAFON(3)
  535. PRINT*,'Y2 :',X(4) ,' DY2 :',DX(4),'ERR :',DELAFON(4)
  536. PRINT*,'Z1 :',X(5) ,' DZ1 :',DX(5),'ERR :',DELAFON(5)
  537. PRINT*,'Z2 :',X(6) ,' DZ2 :',DX(6),'ERR :',DELAFON(6)
  538. PRINT*,'H1',H1,'H2',H2,ht,ht2
  539. PRINT*
  540. write(*,200)jacob
  541. END IF
  542.  
  543. END DO
  544. IF (ERREUR .GT. TEST ) THEN
  545. CONV=.FALSE.
  546. ELSE
  547. CONV=.TRUE.
  548. SIG1= X(1)
  549. SIG2= X(2)
  550. Y1 = X(3)
  551. Y2 = X(4)
  552. Z1 = X(5)
  553. Z2 = X(6)
  554. IF (SIG1 .GE. 0.D0) THEN
  555. JSIGN(1)=1
  556. ELSE
  557. JSIGN(1)=0
  558. END IF
  559. IF (SIG2.GE. 0.D0)THEN
  560. JSIGN(2)=1
  561. ELSE
  562. JSIGN(2)=0
  563. END IF
  564. JSIGN(3)=1
  565. END IF
  566. RETURN
  567. 200 format(6(6(2x,e11.4),/))
  568. END
  569.  
  570.  
  571.  
  572.  
  573.  
  574.  
  575.  
  576.  
  577.  
  578.  
  579.  
  580.  

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