Télécharger modjfd.eso

Retour à la liste

Numérotation des lignes :

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

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