Télécharger cmodcl.eso

Retour à la liste

Numérotation des lignes :

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

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