Télécharger modclb.eso

Retour à la liste

Numérotation des lignes :

  1. C MODCLB SOURCE CB215821 16/04/21 21:17:46 8920
  2. C
  3. C
  4. C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  5. C % %
  6. C % INTEGRATION DU MODCLB UNILATERAL AVEC FERMETURE LINEAIRE %
  7. C % DE FISSURES EN CONTRAINTES PLANES %
  8. C % (MODCLB NON ASSOCIE) %
  9. C % A L"AIDE DE LA METHODE DE NEWTON %
  10. C % %
  11. C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  12. C
  13. SUBROUTINE MODCLB(WRK0,NMATT,EPSP1,EPSP2,YLIM1,YLIM2,Z1,Z2,Y1,Y2,
  14. 1 JSIGN,SIG1,SIG2,CONV)
  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
  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 DELSIG
  42. INTEGER NDIM
  43. PARAMETER (NDIM=6)
  44. REAL*8 JACOB(NDIM,NDIM),J15,J16,J25,J26
  45. REAL*8 DELAFO(NDIM)
  46. REAL*8 X(NDIM),DX(NDIM),X0(NDIM),XX,XIG1,XIG2
  47. REAL*8 ERREUR,TEST
  48. REAL*8 H1,H2,HT,HT2,trf
  49. REAL*8 ALFA1,ALFA2
  50. REAL*8 SIGY
  51. REAL*8 UN,DEUX
  52. REAL*8 ZMAX
  53. INTEGER ITER ,I,DIMPB
  54. LOGICAL LDOM1,LDOM2
  55. SEGMENT WRK0
  56. REAL*8 XMAT(NMATT)
  57. ENDSEGMENT
  58. UN=1.D0
  59. DEUX=2.D0
  60. ZMAX=1.e+04
  61. E = XMAT(1)
  62. NU = XMAT(2)
  63. Y01 = XMAT(5)
  64. Y02 = XMAT(6)
  65. A1 = XMAT(7)
  66. A2 = XMAT(8)
  67. B1 = XMAT(9)
  68. B2 = XMAT(10)
  69. BETA1 = XMAT(11)
  70. BETA2 = XMAT(12)
  71. SIGF = XMAT(13)
  72. SIGY =SQRT(2.D0*E*Y01)
  73. * DELSIG=SIGY*1.D-4
  74. DELSIG=0.d0
  75. * print*,'E,NU,Y01,Y02,A1,A2,B1,B2,BETA1,BETA2'
  76. * print*,E,NU,Y01,Y02,A1,A2,B1,B2,BETA1,BETA2
  77. TEST=1.D-12
  78. ERREUR=1.D0
  79. ITER=0
  80. X(1)=SIG1
  81. X(2)=SIG2
  82. X(3)=Y1
  83. X(4)=Y2
  84. X(5)=Z1
  85. X(6)=Z2
  86. DO I=1, NDIM
  87. X0(I)=X(I)
  88. END DO
  89. C
  90. C
  91. C
  92. JACOB(1,1) = -UN
  93. JACOB(1,2) = 0.D0
  94. JACOB(1,3) = 0.D0
  95. JACOB(1,4) = 0.D0
  96. JACOB(2,1) = 0.D0
  97. JACOB(2,2) = -UN
  98. JACOB(2,3) = 0.D0
  99. JACOB(2,4) = 0.D0
  100. JACOB(3,3) = -UN
  101. JACOB(3,4) = 0.D0
  102. JACOB(3,6) = 0.D0
  103. JACOB(4,3) = 0.D0
  104. JACOB(4,4) = -UN
  105. JACOB(4,5) = 0.D0
  106. JACOB(5,1) = 0.D0
  107. JACOB(5,2) = 0.D0
  108. JACOB(5,4) = 0.D0
  109. JACOB(5,5) = -UN
  110. JACOB(5,6) = 0.D0
  111. JACOB(6,1) = 0.D0
  112. JACOB(6,2) = 0.D0
  113. JACOB(6,3) = 0.D0
  114. JACOB(6,5) = 0.D0
  115. JACOB(6,6) = -UN
  116. alfa1=UN
  117. alfa2=UN
  118.  
  119. DO WHILE (ERREUR .GT. TEST .AND. ITER .LE. 50)
  120. ITER=ITER+1
  121. IF (X(1) .GE. DELSIG) THEN
  122. H1=UN
  123. ELSE IF (X(1) .LE. -DELSIG) THEN
  124. H1=0.D0
  125. ELSE
  126. H1=(X(1)+DELSIG)/(DEUX*DELSIG)
  127. END IF
  128. IF (X(2) .GE. DELSIG) THEN
  129. H2=UN
  130. ELSE IF (X(2) .LE. -DELSIG) THEN
  131. H2=0.D0
  132. ELSE
  133. H2=(X(2)+DELSIG)/(DEUX*DELSIG)
  134. END IF
  135. c
  136. c smooth transition on ht ht2
  137. c
  138. trf = x(1) + x(2)
  139. IF(trf.GE.delsig)THEN
  140. HT=UN
  141. HT2=0.D0
  142. else if(trf.ge.(-delsig)) then
  143. ht = (trf+delsig)/DEUX/delsig
  144. ht2 = UN-ht
  145. else if(trf.ge.(-sigf+delsig)) then
  146. ht = 0.D+00
  147. ht2 = UN
  148. else if(trf.ge.(-sigf-delsig)) then
  149. ht = 0.D+00
  150. c ht2 =UN- (trf-(sigf-delsig))/(DEUX*sigf)
  151. ht2 = (trf+sigf+delsig)/DEUX/delsig
  152. else
  153. ht = 0.D+00
  154. ht2 = 0.D+00
  155. endif
  156. LDOM1= (X(3) .GT. YLIM1 .AND. X(3).GT. Y01
  157. 1 .AND. X(5) .LT. ZMAX)
  158. LDOM2= (X(4) .GT. YLIM2 .AND. X(4).GT. Y02
  159. 1 .AND. X(6) .LT. ZMAX)
  160. * print*,'y1',x(3),'ylim1',ylim1
  161. * print*,'ldom1',ldom1,'ldom2',ldom2
  162. C =====> PROBLEME ELASTIQUE
  163. C
  164. C
  165. C CALCUL DU JACOBIEN
  166. C
  167. C
  168. JACOB(3,1) = H1*X(1)*X(5)*X(5)/E
  169. C
  170. JACOB(3,2) = H2*X(2)*X(5)*X(5)/E
  171. C
  172. JACOB(4,1) = (UN-H1)*X(1)*X(6)*X(6)/E
  173. C
  174. JACOB(4,2) = (UN-H2)*X(2)*X(6)*X(6)/E
  175. C
  176. C CALCUL DU SECOND MEMBRE
  177. C
  178. XIG1 =((((X(5)-UN)*BETA1*HT2+(X(5)-UN)*BETA1*HT-E*EPSP2+(X(6)-UN)
  179. 1 *BETA2)*NU+BETA1*(X(6)*(H2+X(5)*(UN-H2)-UN)+X(5)**2*H2-X(5)*H2)
  180. 2 *HT2+BETA1*(X(6)*(H2+X(5)*(UN-H2)-UN)+X(5)**2*H2-X(5)*H2)*HT+
  181. 3 BETA2*(X(6)*(X(5)*H2+H2-UN)-X(5)*H2+X(6)**2*(UN-H2))+E*EPSP1*
  182. 4 (X(6)*(H2-UN)-X(5)*H2))*SIGF+(BETA1*E*(X(5)*EPSP2-EPSP2)+
  183. 5 (UN-X(5))*BETA1*E*EPSP1)*HT2)/((NU**2+H1*(X(5)*X(6)*(2*H2-UN)-
  184. 6 X(5)**2*H2+X(6)**2*(UN-H2))-X(5)*X(6)*H2+X(6)**2*(H2-UN))*SIGF+
  185. 7 (DEUX-DEUX*X(5))*BETA1*HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*
  186. 8 (-H2+X(5)*(H2-DEUX)+DEUX))+((X(5)-UN)*X(6)-X(5)**2+X(5))
  187. 9 *BETA1*H1)*HT2)
  188. C
  189. XIG2 =((((X(5)-UN)*BETA1*HT2+(X(5)-UN)*BETA1*HT-E*EPSP1+(X(6)-UN)
  190. 1 *BETA2)*NU+(((UN-X(5))*X(6)+X(5)**2-X(5))*BETA1*H1+(X(5)-UN)*
  191. 2 X(6)*BETA1)*HT2+(((UN-X(5))*X(6)+X(5)**2-X(5))*BETA1*H1+
  192. 3 (X(5)-UN)*X(6)*BETA1)*HT+(E*(X(6)*EPSP2-X(5)*EPSP2)+(-X(6)**2+
  193. 4 (X(5)+UN)*X(6)-X(5))*BETA2)*H1-X(6)*E*EPSP2+(X(6)**2-X(6))*
  194. 5 BETA2)*SIGF+(BETA1*E*(EPSP2-X(5)*EPSP2)+(X(5)-UN)*BETA1*E*EPSP1)
  195. 6 *HT2)/((NU**2+H1*(X(5)*X(6)*(DEUX*H2-UN)-X(5)**2*H2+X(6)**2*
  196. 7 (UN-H2))-X(5)*X(6)*H2+X(6)**2*(H2-UN))*SIGF+(DEUX-DEUX*X(5))*
  197. 8 BETA1*HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*
  198. 9 (H2-DEUX)+DEUX))+((X(5)-UN)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)
  199.  
  200. DELAFO(1) =X(1)-XIG1
  201. C
  202. DELAFO(2) =X(2)-XIG2
  203. C
  204. XX= (H1*X(1)*X(1)+H2*X(2)*X(2))*X(5)*X(5)/(DEUX*E)
  205. DELAFO(3) = X(3)-XX
  206. C
  207. XX=((UN-H1)*X(1)*X(1)+(UN-H2)*X(2)*X(2) )
  208. 1 *X(6)*X(6)/(DEUX*E)
  209. DELAFO(4) = X(4)-XX
  210. C
  211. IF (LDOM1) THEN
  212. C
  213. C =====> ENDOMMAGEMENT EN TRACTION
  214. C
  215. J15 = (((BETA1*HT2+BETA1*HT)*NU+BETA1*(DEUX*X(5)*H2-H2+X(6)*(UN-H2
  216. & ))*HT2+BETA1*(DEUX*X(5)*H2-H2+X(6)*(UN-H2))*HT+BETA2*(X(6)*H2-H
  217. & 2)-E*EPSP1*H2)*SIGF+(BETA1*E*EPSP2-BETA1*E*EPSP1)*HT2)/((NU**2+
  218. & H1*(X(5)*X(6)*(DEUX*H2-UN)-X(5)**2*H2+X(6)**2*(UN-H2))-X(5)*X(6
  219. & )*H2+X(6)**2*(H2-UN))*SIGF+(DEUX-DEUX*X(5))*BETA1*HT2*NU+(BETA1
  220. & *(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-DEUX)+DEUX))+((X(5)-UN
  221. & )*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)-((H1*(X(6)*(DEUX*H2-UN)-DEU
  222. & X*X(5)*H2)-X(6)*H2)*SIGF-DEUX*BETA1*HT2*NU+(BETA1*(-DEUX*X(5)*H
  223. & 2+H2+X(6)*(H2-DEUX))+(X(6)-DEUX*X(5)+UN)*BETA1*H1)*HT2)*((((X(5
  224. & )-UN)*BETA1*HT2+(X(5)-UN)*BETA1*HT-E*EPSP2+(X(6)-UN)*BETA2)*NU+
  225. & BETA1*(X(6)*(H2+X(5)*(UN-H2)-UN)+X(5)**2*H2-X(5)*H2)*HT2+BETA1*
  226. & (X(6)*(H2+X(5)*(UN-H2)-UN)+X(5)**2*H2-X(5)*H2)*HT+BETA2*(X(6)*(
  227. & X(5)*H2+H2-UN)-X(5)*H2+X(6)**2*(UN-H2))+E*EPSP1*(X(6)*(H2-UN)-X
  228. & (5)*H2))*SIGF+(BETA1*E*(X(5)*EPSP2-EPSP2)+(UN-X(5))*BETA1*E*EPS
  229. & P1)*HT2)/((NU**2+H1*(X(5)*X(6)*(DEUX*H2-UN)-X(5)**2*H2+X(6)**2*
  230. & (UN-H2))-X(5)*X(6)*H2+X(6)**2*(H2-UN))*SIGF+(DEUX-DEUX*X(5))*BE
  231. & TA1*HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-DEUX)
  232. & +DEUX))+((X(5)-UN)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)**2
  233. C
  234. C
  235. J25 = (((BETA1*HT2+BETA1*HT)*NU+((-X(6)+DEUX*X(5)-UN)*BETA1*H1+X(6
  236. & )*BETA1)*HT2+((-X(6)+DEUX*X(5)-UN)*BETA1*H1+X(6)*BETA1)*HT+((X(
  237. & 6)-UN)*BETA2-E*EPSP2)*H1)*SIGF+(BETA1*E*EPSP1-BETA1*E*EPSP2)*HT
  238. & 2)/((NU**2+H1*(X(5)*X(6)*(DEUX*H2-UN)-X(5)**2*H2+X(6)**2*(UN-H2
  239. & ))-X(5)*X(6)*H2+X(6)**2*(H2-UN))*SIGF+(DEUX-DEUX*X(5))*BETA1*HT
  240. & 2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-DEUX)+DEUX)
  241. & )+((X(5)-UN)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)-((H1*(X(6)*(DEUX
  242. & *H2-UN)-DEUX*X(5)*H2)-X(6)*H2)*SIGF-DEUX*BETA1*HT2*NU+(BETA1*(-
  243. & DEUX*X(5)*H2+H2+X(6)*(H2-DEUX))+(X(6)-DEUX*X(5)+UN)*BETA1*H1)*H
  244. & T2)*((((X(5)-UN)*BETA1*HT2+(X(5)-UN)*BETA1*HT-E*EPSP1+(X(6)-UN)
  245. & *BETA2)*NU+(((UN-X(5))*X(6)+X(5)**2-X(5))*BETA1*H1+(X(5)-UN)*X(
  246. & 6)*BETA1)*HT2+(((UN-X(5))*X(6)+X(5)**2-X(5))*BETA1*H1+(X(5)-UN)
  247. & *X(6)*BETA1)*HT+(E*(X(6)*EPSP2-X(5)*EPSP2)+(-X(6)**2+(X(5)+UN)*
  248. & X(6)-X(5))*BETA2)*H1-X(6)*E*EPSP2+(X(6)**2-X(6))*BETA2)*SIGF+(B
  249. & ETA1*E*(EPSP2-X(5)*EPSP2)+(X(5)-UN)*BETA1*E*EPSP1)*HT2)/((NU**2
  250. & +H1*(X(5)*X(6)*(DEUX*H2-UN)-X(5)**2*H2+X(6)**2*(UN-H2))-X(5)*X(
  251. & 6)*H2+X(6)**2*(H2-UN))*SIGF+(DEUX-DEUX*X(5))*BETA1*HT2*NU+(BETA
  252. & 1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-DEUX)+DEUX))+((X(5)-U
  253. & N)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)**2
  254. C
  255. JACOB(1,5) = J15
  256. JACOB(2,5) = J25
  257. JACOB(3,5) = (H1*X(1)*X(1)+H2*X(2)*X(2))*X(5)/E
  258. JACOB(4,5) = 0.D0
  259. JACOB(5,5) = -1.D0
  260. JACOB(5,1) = 0.D0
  261. JACOB(5,2) = 0.D0
  262. JACOB(5,3) = A1*B1*(A1*(X(3)-Y01))**(B1-UN)
  263. JACOB(5,4) = 0.D0
  264. DELAFO(5) = X(5)-((A1*(X(3)-Y01))**B1+UN)
  265. IF (LDOM2) THEN
  266. C
  267. C =====> ENDOMMAGEMENT EN COMPRESSION
  268. C
  269. J16 = (BETA2*NU+BETA1*(H2+X(5)*(1-H2)-1)*HT2+BETA1*(H2+X(5)*(1-H2)
  270. 1 -1)*HT+BETA2*(X(5)*H2+H2+2*X(6)*(1-H2)-1)+E*EPSP1*(H2-1))*SIGF/
  271. 2 ((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X(5)**2*H2+X(6)**2*(1-H2))-X(5)*
  272. 3 X(6)*H2+X(6)**2*(H2-1))*SIGF+(2-2*X(5))*BETA1*HT2*NU+(BETA1*(-X
  273. 4 (5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-2)+2))+((X(5)-1)*X(6)-X(5)
  274. 5 **2+X(5))*BETA1*H1)*HT2)-((H1*(X(5)*(2*H2-1)+2*X(6)*(1-H2))-X(5
  275. 6 )*H2+2*X(6)*(H2-1))*SIGF+(BETA1*(-H2+X(5)*(H2-2)+2)+(X(5)-1)*BE
  276. 7 TA1*H1)*HT2)*((((X(5)-1)*BETA1*HT2+(X(5)-1)*BETA1*HT-E*EPSP2+(X
  277. 8 (6)-1)*BETA2)*NU+BETA1*(X(6)*(H2+X(5)*(1-H2)-1)+X(5)**2*H2-X(5)
  278. 9 *H2)*HT2+BETA1*(X(6)*(H2+X(5)*(1-H2)-1)+X(5)**2*H2-X(5)*H2)*HT+
  279. : BETA2*(X(6)*(X(5)*H2+H2-1)-X(5)*H2+X(6)**2*(1-H2))+E*EPSP1*(X(6
  280. ; )*(H2-1)-X(5)*H2))*SIGF+(BETA1*E*(X(5)*EPSP2-EPSP2)+(1-X(5))*BE
  281. < TA1*E*EPSP1)*HT2)/((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X(5)**2*H2+X(6
  282. = )**2*(1-H2))-X(5)*X(6)*H2+X(6)**2*(H2-1))*SIGF+(2-2*X(5))*BETA1
  283. > *HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-2)+2))+(
  284. ? (X(5)-1)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)**2
  285. C
  286. J26 = (BETA2*NU+((1-X(5))*BETA1*H1+(X(5)-1)*BETA1)*HT2+((1-X(5))*B
  287. 1 ETA1*H1+(X(5)-1)*BETA1)*HT+(E*EPSP2+(-2*X(6)+X(5)+1)*BETA2)*H1-
  288. 2 E*EPSP2+(2*X(6)-1)*BETA2)*SIGF/((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X
  289. 3 (5)**2*H2+X(6)**2*(1-H2))-X(5)*X(6)*H2+X(6)**2*(H2-1))*SIGF+(2-
  290. 4 2*X(5))*BETA1*HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)
  291. 5 *(H2-2)+2))+((X(5)-1)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)-((H1*(X
  292. 6 (5)*(2*H2-1)+2*X(6)*(1-H2))-X(5)*H2+2*X(6)*(H2-1))*SIGF+(BETA1*
  293. 7 (-H2+X(5)*(H2-2)+2)+(X(5)-1)*BETA1*H1)*HT2)*((((X(5)-1)*BETA1*H
  294. 8 T2+(X(5)-1)*BETA1*HT-E*EPSP1+(X(6)-1)*BETA2)*NU+(((1-X(5))*X(6)
  295. 9 +X(5)**2-X(5))*BETA1*H1+(X(5)-1)*X(6)*BETA1)*HT2+(((1-X(5))*X(6
  296. : )+X(5)**2-X(5))*BETA1*H1+(X(5)-1)*X(6)*BETA1)*HT+(E*(X(6)*EPSP2
  297. ; -X(5)*EPSP2)+(-X(6)**2+(X(5)+1)*X(6)-X(5))*BETA2)*H1-X(6)*E*EPS
  298. < P2+(X(6)**2-X(6))*BETA2)*SIGF+(BETA1*E*(EPSP2-X(5)*EPSP2)+(X(5)
  299. = -1)*BETA1*E*EPSP1)*HT2)/((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X(5)**2*
  300. > H2+X(6)**2*(1-H2))-X(5)*X(6)*H2+X(6)**2*(H2-1))*SIGF+(2-2*X(5))
  301. ? *BETA1*HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-2)
  302. @ +2))+((X(5)-1)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)**2
  303. C
  304. C
  305. JACOB(1,6) = J16
  306. JACOB(2,6) = J26
  307. JACOB(3,6) = 0.D0
  308. JACOB(4,6) = ((1.D0-H1)*X(1)*X(1)+(1.D0-H2)
  309. & *X(2)*X(2))*X(6)/E
  310. JACOB(5,6) = 0.D0
  311. JACOB(6,6) = -1.D0
  312. JACOB(6,1) = 0.D0
  313. JACOB(6,2) = 0.D0
  314. JACOB(6,3) = 0.D0
  315. * JACOB(6,4) = (A2*B2*(X(4)-Y02)**(B2-1.0))
  316. JACOB(6,4) = A2*B2*(A2*(X(4)-Y02))**(B2-1.0)
  317. JACOB(6,5) = 0.D0
  318. * DELAFO(6) = X(6)-(A2*(X(4)-Y02)**B2+1.0)
  319. DELAFO(6) = X(6)-((A2*(X(4)-Y02))**B2+1.0)
  320. DIMPB=6
  321. ELSE
  322. DIMPB=5
  323. END IF
  324. ELSE
  325. C
  326. C =====> PAS D'ENDOMMAGEMENT EN TRACTION
  327. C
  328. IF (LDOM2) THEN
  329. C
  330. C =====> ENDOMMAGEMENT EN COMPRESSION
  331. C
  332. J16 = (BETA2*NU+BETA1*(H2+X(5)*(1-H2)-1)*HT2+BETA1*(H2+X(5)*(1-H2)
  333. 1 -1)*HT+BETA2*(X(5)*H2+H2+2*X(6)*(1-H2)-1)+E*EPSP1*(H2-1))*SIGF/
  334. 2 ((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X(5)**2*H2+X(6)**2*(1-H2))-X(5)*
  335. 3 X(6)*H2+X(6)**2*(H2-1))*SIGF+(2-2*X(5))*BETA1*HT2*NU+(BETA1*(-X
  336. 4 (5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-2)+2))+((X(5)-1)*X(6)-X(5)
  337. 5 **2+X(5))*BETA1*H1)*HT2)-((H1*(X(5)*(2*H2-1)+2*X(6)*(1-H2))-X(5
  338. 6 )*H2+2*X(6)*(H2-1))*SIGF+(BETA1*(-H2+X(5)*(H2-2)+2)+(X(5)-1)*BE
  339. 7 TA1*H1)*HT2)*((((X(5)-1)*BETA1*HT2+(X(5)-1)*BETA1*HT-E*EPSP2+(X
  340. 8 (6)-1)*BETA2)*NU+BETA1*(X(6)*(H2+X(5)*(1-H2)-1)+X(5)**2*H2-X(5)
  341. 9 *H2)*HT2+BETA1*(X(6)*(H2+X(5)*(1-H2)-1)+X(5)**2*H2-X(5)*H2)*HT+
  342. : BETA2*(X(6)*(X(5)*H2+H2-1)-X(5)*H2+X(6)**2*(1-H2))+E*EPSP1*(X(6
  343. ; )*(H2-1)-X(5)*H2))*SIGF+(BETA1*E*(X(5)*EPSP2-EPSP2)+(1-X(5))*BE
  344. < TA1*E*EPSP1)*HT2)/((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X(5)**2*H2+X(6
  345. = )**2*(1-H2))-X(5)*X(6)*H2+X(6)**2*(H2-1))*SIGF+(2-2*X(5))*BETA1
  346. > *HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-2)+2))+(
  347. ? (X(5)-1)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)**2
  348. C
  349. J26 = (BETA2*NU+((1-X(5))*BETA1*H1+(X(5)-1)*BETA1)*HT2+((1-X(5))*B
  350. 1 ETA1*H1+(X(5)-1)*BETA1)*HT+(E*EPSP2+(-2*X(6)+X(5)+1)*BETA2)*H1-
  351. 2 E*EPSP2+(2*X(6)-1)*BETA2)*SIGF/((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X
  352. 3 (5)**2*H2+X(6)**2*(1-H2))-X(5)*X(6)*H2+X(6)**2*(H2-1))*SIGF+(2-
  353. 4 2*X(5))*BETA1*HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)
  354. 5 *(H2-2)+2))+((X(5)-1)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)-((H1*(X
  355. 6 (5)*(2*H2-1)+2*X(6)*(1-H2))-X(5)*H2+2*X(6)*(H2-1))*SIGF+(BETA1*
  356. 7 (-H2+X(5)*(H2-2)+2)+(X(5)-1)*BETA1*H1)*HT2)*((((X(5)-1)*BETA1*H
  357. 8 T2+(X(5)-1)*BETA1*HT-E*EPSP1+(X(6)-1)*BETA2)*NU+(((1-X(5))*X(6)
  358. 9 +X(5)**2-X(5))*BETA1*H1+(X(5)-1)*X(6)*BETA1)*HT2+(((1-X(5))*X(6
  359. : )+X(5)**2-X(5))*BETA1*H1+(X(5)-1)*X(6)*BETA1)*HT+(E*(X(6)*EPSP2
  360. ; -X(5)*EPSP2)+(-X(6)**2+(X(5)+1)*X(6)-X(5))*BETA2)*H1-X(6)*E*EPS
  361. < P2+(X(6)**2-X(6))*BETA2)*SIGF+(BETA1*E*(EPSP2-X(5)*EPSP2)+(X(5)
  362. = -1)*BETA1*E*EPSP1)*HT2)/((NU**2+H1*(X(5)*X(6)*(2*H2-1)-X(5)**2*
  363. > H2+X(6)**2*(1-H2))-X(5)*X(6)*H2+X(6)**2*(H2-1))*SIGF+(2-2*X(5))
  364. ? *BETA1*HT2*NU+(BETA1*(-X(5)**2*H2+X(5)*H2+X(6)*(-H2+X(5)*(H2-2)
  365. @ +2))+((X(5)-1)*X(6)-X(5)**2+X(5))*BETA1*H1)*HT2)**2
  366. C
  367. C
  368.  
  369. JACOB(1,5) = J16
  370. JACOB(2,5) = J26
  371. JACOB(3,5) = 0.D0
  372. JACOB(4,5) = ((1.D0-H1)*X(1)*X(1)+(1.D0-H2)
  373. & *X(2)*X(2))*X(6)/E
  374. JACOB(5,5) = -1.D0
  375. JACOB(5,1) = 0.D0
  376. JACOB(5,2) = 0.D0
  377. JACOB(5,3) = 0.D0
  378. * JACOB(5,4) = (A2*B2*(X(4)-Y02)**(B2-1.0))
  379. JACOB(5,4) = A2*B2*(A2*(X(4)-Y02))**(B2-1.0)
  380. JACOB(5,5) = -1.D0
  381. * DELAFO(5) = X(6)-(A2*(X(4)-Y02)**B2+1.0)
  382. DELAFO(5) = X(6)-((A2*(X(4)-Y02))**B2+1.0)
  383. DIMPB=5
  384. ELSE
  385. DIMPB=4
  386. END IF
  387. END IF
  388. C
  389. C RESOLUTION
  390. C
  391. CALL GAUSCL(JACOB,DX,DELAFO,NDIM,DIMPB,0)
  392. C
  393. C CALCUL DE L"ERREUR ET DE LA NOUVELLE SOLUTION
  394. C
  395. ERREUR=0.D0
  396. DO I=1,2
  397. ERREUR=ERREUR+DX(I)*DX(I)/(SIGY*SIGY)
  398. X(I)=X(I)+DX(I)
  399. END DO
  400. X(3)=X(3)+DX(3)
  401. ERREUR=ERREUR+DX(3)*DX(3)/(Y01*Y01)
  402. X(4)=X(4)+DX(4)
  403. ERREUR=ERREUR+DX(4)*DX(4)/(Y02*Y02)
  404. IF (LDOM1) THEN
  405. IF ((X(5)+DX(5)).LT.Z1) THEN
  406. DX(5)=Z1-X(5)
  407. END IF
  408. IF ( X(5) .GT. ZMAX)THEN
  409. DX(5)=ZMAX + 0.1D0 - X(5)
  410. END IF
  411. X(5)=X(5)+DX(5)
  412. ERREUR=ERREUR+DX(5)*DX(5)/(X(5)*X(5))
  413. IF (LDOM2) THEN
  414. IF ((X(6)+DX(6)).LT.Z2) THEN
  415. DX(6)=Z2-X(6)
  416. END IF
  417. IF ( X(6) .GT. ZMAX)THEN
  418. DX(6)=ZMAX + 0.1D0 - X(6)
  419. END IF
  420. X(6)=X(6)+DX(6)
  421. ERREUR=ERREUR+DX(6)*DX(6)/(X(6)*X(6))
  422. END IF
  423. ELSE IF (LDOM2) THEN
  424. IF ((X(6)+DX(5)).LT.Z2) THEN
  425. DX(5)=Z2-X(6)
  426. END IF
  427. IF ( X(6) .GT. ZMAX)THEN
  428. DX(5)=ZMAX + 0.1D0 - X(6)
  429. END IF
  430. X(6)=X(6)+DX(5)
  431. ERREUR=ERREUR+DX(5)*DX(5)/(X(6)*X(6))
  432. END IF
  433. IF (ERREUR .GT. 1.D+100)THEN
  434. CONV=.FALSE.
  435. print*,'erreur dans modele',erreur
  436. print*,'iteration num',iter
  437. print*,'delafo',delafo
  438. RETURN
  439. END IF
  440. C
  441. C QUELQUES IMPRESSIONS EN CAS D'ERREUR
  442. C
  443. if (.not. conv )then
  444. PRINT*,'ITERATION ',ITER
  445. print*,'LDOM1 :',LDOM1,' LDOM2 :',LDOM2
  446. PRINT*,'EPS1= ',EPSP1,' EPS2= ',EPSP2
  447. PRINT*,'SIG1:',X(1) ,' DSIG1:',DX(1),'ERR :',DELAFO(1)
  448. PRINT*,'SIG2:',X(2) ,' DSIG2:',DX(2),'ERR :',DELAFO(2)
  449. PRINT*,'Y1 :',X(3) ,' DY1 :',DX(3),'ERR :',DELAFO(3)
  450. PRINT*,'Y2 :',X(4) ,' DY2 :',DX(4),'ERR :',DELAFO(4)
  451. PRINT*,'Z1 :',X(5) ,' DZ1 :',DX(5),'ERR :',DELAFO(5)
  452. PRINT*,'Z2 :',X(6) ,' DZ2 :',DX(6),'ERR :',DELAFO(6)
  453. PRINT*,'H1',H1,'H2',H2,ht,ht2
  454. PRINT*
  455. write(*,200)jacob
  456. 200 format(6(6(2x,e11.4),/))
  457. END IF
  458.  
  459. END DO
  460. IF (ERREUR .GT. TEST ) THEN
  461. CONV=.FALSE.
  462. ELSE
  463. CONV=.TRUE.
  464. SIG1= X(1)
  465. SIG2= X(2)
  466. Y1 = X(3)
  467. Y2 = X(4)
  468. Z1 = X(5)
  469. Z2 = X(6)
  470. IF (SIG1 .GE. 0.D0) THEN
  471. JSIGN(1)=1
  472. ELSE
  473. JSIGN(1)=0
  474. END IF
  475. IF (SIG2.GE. 0.D0)THEN
  476. JSIGN(2)=1
  477. ELSE
  478. JSIGN(2)=0
  479. END IF
  480. JSIGN(3)=1
  481. END IF
  482. RETURN
  483. END
  484.  
  485.  
  486.  
  487.  

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