Télécharger modjf2.eso

Retour à la liste

Numérotation des lignes :

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

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