Télécharger mshear.eso

Retour à la liste

Numérotation des lignes :

  1. C MSHEAR SOURCE PV 11/03/07 21:17:33 6885
  2. SUBROUTINE MSHEAR(DDEP,FOR0,FORF,VAR0,VARF,NVAR,DDINL,
  3. & DTRANP,DTRANM,BETA,NPELA,TRFA,DOCP,DOCM,EXPN,
  4. & CURFP,NCURFP,CURKP,NCURKP,CURLP,NCURLP,
  5. & CURFM,NCURFM,CURKM,NCURKM,CURLM,NCURLM,KERRE)
  6. IMPLICIT INTEGER(I-N)
  7. IMPLICIT REAL*8(A-H,O-Z)
  8. C
  9. C======================================================================
  10. C MUR_SHEAR, Elisa, Armelle, Alessandra et Pierre, ELSA/ISPRA
  11. C 01/2002 --> ../2004 ...
  12. C======================================================================
  13. C
  14. C MODELE GLOBAL DE MUR EN CISAILLEMENT
  15. C (Sur des elements de poutre TIMO - Effort tranchant/Cisail.)
  16. C
  17. C======================================================================
  18. C=======================================================================
  19. C
  20. C LISTE D'ECHANGE
  21. C ---------------
  22. C
  23. C DDEP = Incrément de déplacement
  24. C FOR0 = Effort initial
  25. C FORF = Effort final
  26. C VAR0 = Variables internes initiales
  27. C VARF = Variables internes finales
  28. C NVAR = Nb de variable interne
  29. C DDINL = Increment de deformation inelastique
  30. C
  31. C DTRANP = zone de transition (sens positif)
  32. C DTRANM = (sens négatif)
  33. C
  34. C BETA = Coefficient BETA (equivalence nrj-deplacement)
  35. C NPELA = Nb de pts par unite de zone elastique pour les
  36. C segments non lineaires
  37. C+2003
  38. C TRFA = Coefficient d'entrainement des origine (>0,<1)
  39. C DOCP = Deplacement d'ouverture des fissure + (>0)
  40. C DOCM = Deplacement d'ouverture des fissure - (<0)
  41. C+2003
  42. C+2004
  43. C EXPN = Exposant n de la fonction de cumul du domage
  44. C+2004
  45. C
  46. C CURFP,NCURFP = courbe de charge + (x>0,y>0)
  47. C CURKP,NCURKP = courbe de raideur + (x>0,y>0)
  48. C CURLP,NCURLP = courbe de domaine elastique + (x>0,y>0)
  49. C CURFM,NCURFM = courbe de charge - (x<0,y<0)
  50. C CURKM,NCURKM = courbe de raideur - (x<0,y>0)
  51. C CURLM,NCURLM = courbe de domaine elastique - (x<0,y>0)
  52. C
  53. C VARIABLES INTERNES (0=initiales, F=finales)
  54. C ------------------
  55. C
  56. C : Numero du cas
  57. C DEP : Deplacement total
  58. C EPLUS : Energie "positive"
  59. C EMOIN : Energie "negative"
  60. C DPLUS : Deplacement de reference +
  61. C DMOIN : Deplacement de reference -
  62. C FCINI : Force de reference pour la surface cinematique
  63. C FCAMP : Taille de la surface cinematique
  64. C KCINE : Raideur de la surface cinematique
  65. C APLUS : Coefficient multiplicateur de l'amplitude +
  66. C AMOIN : Coefficient multiplicateur de l'amplitude -
  67. C+2003
  68. C OPLUS : Position de l'origine + <0
  69. C OMOIN : Position de l'origine - <0
  70. C+2003
  71. C
  72. C
  73. C=======================================================================
  74. C
  75. PARAMETER (XZER=0.D0,UN=1.D0,XPETIT=1.D-7)
  76. DIMENSION CURFP(2,NCURFP),CURKP(2,NCURKP),CURLP(2,NCURLP)
  77. DIMENSION CURFM(2,NCURFM),CURKM(2,NCURKM),CURLM(2,NCURLM)
  78. C
  79. REAL*8 KCINE0,KCINEF
  80. C
  81. DIMENSION VAR0(NVAR),VARF(NVAR)
  82. LOGICAL LSWICH
  83. C
  84. KERRE=0
  85. C
  86. C Lecture variables internes
  87. C
  88. ICAS = nint(VAR0(1))
  89. DEP0 = VAR0(2)
  90. EPLUS0= VAR0(3)
  91. EMOIN0= VAR0(4)
  92. DPLUS0= VAR0(5)
  93. DMOIN0= VAR0(6)
  94. FCINI0= VAR0(7)
  95. FCAMP0= VAR0(8)
  96. KCINE0= VAR0(9)
  97. APLUS0= VAR0(10)
  98. AMOIN0= VAR0(11)
  99. C+2003
  100. OPLUS0= VAR0(12)
  101. OMOIN0= VAR0(13)
  102. C+2003
  103. C+2004
  104. E= VAR0(12)
  105. OMOIN0= VAR0(13)
  106. C+2004
  107. C
  108. C Initialisation des variables finales
  109. C
  110. DDINL=XZER
  111. DEPF =DEP0+DDEP
  112. DO IE1=1,NVAR
  113. VARF(IE1)=VAR0(IE1)
  114. ENDDO
  115. VARF(2)=DEPF
  116. FORF=FOR0
  117. *
  118. * WRITE(6,*)ICAS,DEP0,DEPF
  119. *
  120. C
  121. C Cas de l'increment nul
  122. C
  123. IF(DDEP.EQ.XZER)RETURN
  124. C==================================================================
  125. C
  126. C DRIVER
  127. C
  128. C==================================================================
  129. GOTO(1,2,3,4,5,6,7),ICAS+1
  130. C==========================================================
  131. C A: CAS DU MODELE DANS SON ETAT ELASTIQUE ICAS= 0
  132. C==========================================================
  133. 1 CONTINUE
  134. C-->Condition de branchement sur B0 ou C0
  135. IF(DEPF.GT.CURFP(1,2))THEN
  136. ICAS=1
  137. DEP0 = CURFP(1,2)
  138. FOR0 = CURFP(2,2)
  139. KCINE0 = FOR0/DEP0
  140. DMOIN0 = CURFM(1,2)
  141. VARF(6)= DMOIN0
  142. GOTO 2
  143. ELSEIF(DEPF.LT.CURFM(1,2))THEN
  144. ICAS=4
  145. DEP0 = CURFM(1,2)
  146. FOR0 = CURFM(2,2)
  147. KCINE0 = FOR0/DEP0
  148. DPLUS0 = CURFP(1,2)
  149. VARF(5)= DPLUS0
  150. GOTO 5
  151. ENDIF
  152. C-->Regime normal: On calcul l'etat elastique
  153. IF(DEPF.GE.0)THEN
  154. KCINEF=CURFP(2,2)/CURFP(1,2)
  155. FORF=KCINEF*DEPF
  156. ELSE
  157. KCINEF=CURFM(2,2)/CURFM(1,2)
  158. FORF=KCINEF*DEPF
  159. ENDIF
  160. RETURN
  161. C==========================================================
  162. C B0: CAS DU CHARGEMENT MONOTONE DIRECTION + ICAS= 1
  163. C==========================================================
  164. 2 CONTINUE
  165. C-->Condition de branchement sur C1
  166. IF(DDEP.LT.XZER)THEN
  167. ICAS=2
  168. GOTO 3
  169. ENDIF
  170. C-->Regime normal: On procede sur la courbe de charge +
  171. C 2003DPLUSF=DEPF
  172. DPLUSF=DEPF-OPLUS0
  173. C 2003CALL YOFXCU(DEPF ,CURFP,NCURFP, FORF,DYSDX,KERRE)
  174. CALL YOFXCU(DPLUSF,CURFP,NCURFP, FORF,DYSDX,KERRE)
  175. FORF=FORF*(UN-APLUS0)
  176. VARF(1)=ICAS
  177. VARF(5)=DPLUSF
  178. C-->On calcule et on sauve les caracteristiques de la surface
  179. C Cinematique en tenant compte de la zone de transition
  180. C 2003CALL MSHEA1(DEPF,ICAS,FORF,DPLUSF,DMOIN0, DTRANP,DTRANM,
  181. CALL MSHEA1(DPLUSF,ICAS,FORF,DPLUSF,DMOIN0, DTRANP,DTRANM,
  182. > APLUS0,CURFP,NCURFP,CURKP,NCURKP,CURLP,NCURLP,
  183. > AMOIN0,CURFM,NCURFM,CURKM,NCURKM,CURLM,NCURLM,
  184. > FCINIF,FCAMPF,KCINEF)
  185. C2003> FCINIF,FCAMPF,KCINEF,ALPHP)
  186. VARF(7) =FCINIF
  187. VARF(8) =FCAMPF
  188. VARF(9) =KCINEF
  189. C-->Increment de d.p.
  190. DDINL=DDINL+DEPF-FORF/KCINEF-(DEP0-FOR0/KCINE0)
  191. C+2003
  192. IF(DPLUSF.GE.DOCP)THEN
  193. IF(OMOIN0.EQ.XZER)THEN
  194. OMOINF=TRFA*(DPLUSF-DOCP)
  195. ELSE
  196. OMOINF=OMOIN0+TRFA*(DPLUSF-DPLUS0)
  197. ENDIF
  198. VARF(13)=OMOINF
  199. ENDIF
  200. C+2003
  201. RETURN
  202. C==========================================================
  203. C C1: BRANCHE DE UNLOADING ELASTIQUE - ICAS= 2
  204. C==========================================================
  205. 3 CONTINUE
  206. C-->Condition de branchement sur B0, B2 ou C2
  207. FORF=FOR0+KCINE0*DDEP
  208. IF(FORF.GT.FCINI0)THEN
  209. C -->On sort par le haut (B0 ou B2)
  210. DEP0=DEP0+(FCINI0-FOR0)/KCINE0
  211. DDEP=DEPF-DEP0
  212. FOR0=FCINI0
  213. C -->On retourne sur B0... si on vient de B0
  214. IF(ABS(DEP0-DPLUS0).LT.XPETIT*DPLUS0)THEN
  215. ICAS=1
  216. GOTO 2
  217. C -->On retourne sur B2... si on ne vient PAS de B0
  218. ELSE
  219. ICAS=6
  220. GOTO 7
  221. ENDIF
  222. ELSEIF(FORF.LT.FCINI0-FCAMP0)THEN
  223. C --> On continue sur C2
  224. DEP0=DEP0+(FCINI0-FCAMP0-FOR0)/KCINE0
  225. DDEP=DEPF-DEP0
  226. FOR0=FCINI0-FCAMP0
  227. C
  228. C+2002 On repointe "tangent" a la courbe (si possible)
  229. CALL YOFXCU(DMOIN0,CURFM,NCURFM, FMOIN0,DYSDX,KERRE)
  230. C 2003 STRIAL=(FMOIN0-FOR0)/(DMOIN0-DEP0)
  231. STRIAL=(FMOIN0-FOR0)/(DMOIN0+OMOIN0-DEP0)
  232. DO IE1=1,NCURFM
  233. IF(DMOIN0.GT.CURFM(1,IE1))THEN
  234. C 2003 SNEXT=(CURFM(2,IE1)-FOR0)/(CURFM(1,IE1)-DEP0)
  235. SNEXT=(CURFM(2,IE1)-FOR0)/(CURFM(1,IE1)+OMOIN0-DEP0)
  236. IF(SNEXT.GT.STRIAL)THEN
  237. DMOIN0=CURFM(1,IE1)
  238. STRIAL=SNEXT
  239. ELSE
  240. GOTO 31
  241. ENDIF
  242. ENDIF
  243. ENDDO
  244. 31 CONTINUE
  245. VARF(6)=DMOIN0
  246. C+2002
  247. ICAS=3
  248. GOTO 4
  249. ENDIF
  250. C-->Regime normal: On decharge elastiquement
  251. VARF(1)=ICAS
  252. RETURN
  253. C==========================================================
  254. C C2: BRANCHE DE UNLOADING ANELASTIQUE - ICAS= 3
  255. C==========================================================
  256. 4 CONTINUE
  257. C-->Condition de branchement sur B1
  258. IF(DDEP.GT.XZER)THEN
  259. ICAS=5
  260. GOTO 6
  261. ENDIF
  262. C-->Valeur de ALPHP
  263. C 2003ALPHP=(DEP0-FOR0/KCINE0-DTRANM)/(DTRANP-DTRANM)
  264. C 2003ALPHP=MAX(XZER,ALPHP)
  265. C 2003ALPHP=MIN(UN ,ALPHP)
  266. C+2003
  267. CC CALL MSHEA2(DEP0,DTRANP+OPLUS0,DTRANM+OMOIN0,
  268. CALL MSHEA2(DEP0,DTRANP+OPLUS0+OMOIN0,DTRANM+OPLUS0+OMOIN0,
  269. > FOR0,FCAMP0,KCINE0,ICAS,
  270. > ALPHP,ALPHM)
  271. C+2003
  272. C-->Determination de la taille du sous increment
  273. C WARNING DDEP et XINCR sont negatifs
  274. XINCR=(CURFM(2,2)-CURFP(2,2))/NPELA
  275. NINCR=INT(DDEP/XINCR)+1
  276. XINCR=DDEP/NINCR
  277. C-->Determination initiale du point "parfait" que l'on vise (DMOIN0,FMOIN0)
  278. CALL YOFXCU(DMOIN0,CURFM,NCURFM, FMOIN0,DYSDX,KERRE)
  279. C-->Loop sur les sous-increments
  280. LSWICH=.FALSE.
  281. DO IE1=1,NINCR
  282. C-->Traitement du branchement possible sur C0
  283. C 2003 IF(DEP0+XINCR.LT.DMOIN0)THEN
  284. C IF(DEP0+XINCR.LT.DMOIN0+OMOIN0)THEN
  285. C 2003 IF(DEP0+XINCR.LT.(UN-XPETIT)*DMOIN0)THEN
  286. IF(DEP0+XINCR.LT.(UN-XPETIT)*(DMOIN0+OMOIN0))THEN
  287. LSWICH=.TRUE.
  288. C 2003 XINCR =DMOIN0-DEP0
  289. XINCR =DMOIN0+OMOIN0-DEP0
  290. ENDIF
  291. C-->On incremente explicitement tout ce qui doit etre incremente
  292. C WARNING:force effective pointee=FMOIN0*(UN-AMOIN0)
  293. C 2003 FORF=FOR0+XINCR/(DMOIN0-DEP0)*(FMOIN0*(UN-AMOIN0)-FOR0)
  294. FORF=FOR0+XINCR/(DMOIN0+OMOIN0-DEP0)*(FMOIN0*(UN-AMOIN0)-FOR0)
  295. DEP0=DEP0+XINCR
  296. DENE=ABS(FCAMP0/2*(XINCR-(FORF-FOR0)/KCINE0))
  297. DDINL=DDINL+XINCR+FOR0/KCINE0
  298. FOR0=FORF
  299. C
  300. EPLUS0=EPLUS0+ ALPHP *DENE
  301. C 2003 EMOIN0=EMOIN0+(UN-ALPHP)*DENE
  302. EMOIN0=EMOIN0+ ALPHM *DENE
  303. * APLUS0=MIN(UN,APLUS0+BETA* ALPHP* DENE)
  304. * AMOIN0=MIN(UN,AMOIN0+BETA*(UN-ALPHP)*DENE)
  305. C 2003 APLUS0=TANH(BETA*EPLUS0)
  306. C 2003 AMOIN0=TANH(BETA*EMOIN0)
  307. CALL MSHEA3(BETA,EXPN, EPLUS0,EMOIN0, APLUS0,AMOIN0)
  308. C-->On calcule les caracteristiques de la surface
  309. C Cinematique en tenant compte de la zone de transition
  310. CALL MSHEA1(DEP0,ICAS,FOR0,DPLUS0,DMOIN0,
  311. > DTRANP,DTRANM,
  312. > APLUS0,CURFP,NCURFP,CURKP,NCURKP,CURLP,NCURLP,
  313. > AMOIN0,CURFM,NCURFM,CURKM,NCURKM,CURLM,NCURLM,
  314. > FCINI0,FCAMP0,KCINE0)
  315. C2003> FCINI0,FCAMP0,KCINE0,ALPHP)
  316. C+2003
  317. CC CALL MSHEA2(DEP0,DTRANP+OPLUS0,DTRANM+OMOIN0,
  318. CALL MSHEA2(DEP0,DTRANP+OPLUS0+OMOIN0,DTRANM+OPLUS0+OMOIN0,
  319. > FOR0,FCAMP0,KCINE0,ICAS,
  320. > ALPHP,ALPHM)
  321. C+2003
  322. DDINL=DDINL-FOR0/KCINE0
  323. C-->Switch sur C0
  324. IF(LSWICH)THEN
  325. ICAS=4
  326. VARF(3)=EPLUS0
  327. VARF(4)=EMOIN0
  328. VARF(10)=APLUS0
  329. VARF(11)=AMOIN0
  330. GOTO 5
  331. ENDIF
  332. C-->Fin normale loop sur les sous-increments
  333. ENDDO
  334. C
  335. FORF=FOR0
  336. C
  337. VARF(1)=ICAS
  338. VARF(3)=EPLUS0
  339. VARF(4)=EMOIN0
  340. VARF(7)=FCINI0
  341. VARF(8)=FCAMP0
  342. VARF(9)=KCINE0
  343. VARF(10)=APLUS0
  344. VARF(11)=AMOIN0
  345. C
  346. RETURN
  347. C==========================================================
  348. C B0: CAS DU CHARGEMENT MONOTONE DIRECTION - ICAS= 4
  349. C==========================================================
  350. 5 CONTINUE
  351. C-->Condition de branchement sur B1
  352. IF(DDEP.GT.XZER)THEN
  353. ICAS=5
  354. GOTO 6
  355. ENDIF
  356. C-->Regime normal: On procede sur la courbe de charge -
  357. C 2003DMOINF=DEPF
  358. DMOINF=DEPF-OMOIN0
  359. C 2003CALL YOFXCU(DEPF ,CURFM,NCURFM, FORF,DYSDX,KERRE)
  360. CALL YOFXCU(DMOINF,CURFM,NCURFM, FORF,DYSDX,KERRE)
  361. FORF=FORF*(UN-AMOIN0)
  362. VARF(1)=ICAS
  363. VARF(6)=DMOINF
  364. C-->On calcule et on sauve les caracteristiques de la surface
  365. C Cinematique en tenant compte de la zone de transition
  366. C 2003CALL MSHEA1(DEPF,ICAS,FORF,DPLUS0,DMOINF, DTRANP,DTRANM,
  367. CALL MSHEA1(DEPF,ICAS,FORF,DPLUS0,DMOINF, DTRANP,DTRANM,
  368. > APLUS0,CURFP,NCURFP,CURKP,NCURKP,CURLP,NCURLP,
  369. > AMOIN0,CURFM,NCURFM,CURKM,NCURKM,CURLM,NCURLM,
  370. > FCINIF,FCAMPF,KCINEF)
  371. C2003> FCINIF,FCAMPF,KCINEF,ALPHP)
  372. VARF(7) =FCINIF
  373. VARF(8) =FCAMPF
  374. VARF(9) =KCINEF
  375. C-->Increment de d.p.
  376. DDINL=DDINL+DEPF-FORF/KCINEF-(DEP0-FOR0/KCINE0)
  377. C+2003
  378. IF(DMOINF.LE.DOCM)THEN
  379. IF(OPLUS0.EQ.XZER)THEN
  380. OPLUSF=TRFA*(DMOINF-DOCM)
  381. ELSE
  382. OPLUSF=OPLUS0+TRFA*(DMOINF-DMOIN0)
  383. ENDIF
  384. VARF(12)=OPLUSF
  385. ENDIF
  386. C+2003
  387. RETURN
  388. C==========================================================
  389. C B1: BRANCHE DE RELOADING ELASTIQUE + ICAS= 5
  390. C==========================================================
  391. 6 CONTINUE
  392. C-->Condition de branchement sur C0, C2 ou B2
  393. FORF=FOR0+KCINE0*DDEP
  394. IF(FORF.LT.FCINI0)THEN
  395. C -->On sort par le bas (C0 ou C2)
  396. DEP0=DEP0+(FCINI0-FOR0)/KCINE0
  397. DDEP=DEPF-DEP0
  398. FOR0=FCINI0
  399. C -->On retourne sur C0... si on vient de C0
  400. IF(ABS(DEP0-DMOIN0).LT.ABS(XPETIT*DMOIN0))THEN
  401. ICAS=4
  402. GOTO 5
  403. C -->On retourne sur C2... si on ne vient PAS de C0
  404. ELSE
  405. ICAS=3
  406. GOTO 4
  407. ENDIF
  408. ELSEIF(FORF.GT.FCINI0+FCAMP0)THEN
  409. C --> On continue sur B2
  410. DEP0=DEP0+(FCINI0+FCAMP0-FOR0)/KCINE0
  411. DDEP=DEPF-DEP0
  412. FOR0=FCINI0+FCAMP0
  413. C
  414. C2002 On repointe "tangent" a la courbe (si possible)
  415. CALL YOFXCU(DPLUS0,CURFP,NCURFP, FPLUS0,DYSDX,KERRE)
  416. C 2003 STRIAL=(FPLUS0-FOR0)/(DPLUS0-DEP0)
  417. STRIAL=(FPLUS0-FOR0)/(DPLUS0+OPLUS0-DEP0)
  418. DO IE1=1,NCURFP
  419. IF(DPLUS0.LT.CURFP(1,IE1))THEN
  420. C 2003 SNEXT=(CURFP(2,IE1)-FOR0)/(CURFP(1,IE1)-DEP0)
  421. SNEXT=(CURFP(2,IE1)-FOR0)/(CURFP(1,IE1)+OPLUS0-DEP0)
  422. IF(SNEXT.GT.STRIAL)THEN
  423. DPLUS0=CURFP(1,IE1)
  424. STRIAL=SNEXT
  425. ELSE
  426. GOTO 61
  427. ENDIF
  428. ENDIF
  429. ENDDO
  430. 61 CONTINUE
  431. VARF(5)=DPLUS0
  432. C2002
  433. ICAS=6
  434. GOTO 7
  435. ENDIF
  436. C-->Regime normal: On decharge elastiquement
  437. VARF(1)=ICAS
  438. RETURN
  439. C==========================================================
  440. C B2: BRANCHE DE UNLOADING ANELASTIQUE - ICAS= 6
  441. C==========================================================
  442. 7 CONTINUE
  443. C-->Condition de branchement sur C1
  444. IF(DDEP.LT.XZER)THEN
  445. ICAS=2
  446. GOTO 3
  447. ENDIF
  448. C-->Valeur de ALPHP
  449. C 2003ALPHP=(DEP0-FOR0/KCINE0-DTRANM)/(DTRANP-DTRANM)
  450. C 2003ALPHP=MAX(XZER,ALPHP)
  451. C 2003ALPHP=MIN(UN ,ALPHP)
  452. C+2003
  453. CC CALL MSHEA2(DEP0,DTRANP+OPLUS0,DTRANM+OMOIN0,
  454. CALL MSHEA2(DEP0,DTRANP+OPLUS0+OMOIN0,DTRANM+OPLUS0+OMOIN0,
  455. > FOR0,FCAMP0,KCINE0,ICAS,
  456. > ALPHP,ALPHM)
  457. C+2003
  458.  
  459. C-->Determination de la taille du sous increment
  460. C WARNING DDEP et XINCR sont negatifs
  461. XINCR=(CURFP(2,2)-CURFM(2,2))/NPELA
  462. NINCR=INT(DDEP/XINCR)+1
  463. XINCR=DDEP/NINCR
  464. C-->Determination initiale du point "parfait" que l'on vise (DPLUS0,FPLUS0)
  465. CALL YOFXCU(DPLUS0,CURFP,NCURFP, FPLUS0,DYSDX,KERRE)
  466. C-->Loop sur les sous-increments
  467. LSWICH=.FALSE.
  468. DO IE1=1,NINCR
  469. C-->Traitement du branchement possible sur B0
  470. C 2003 IF(DEP0+XINCR.GT.DPLUS0)THEN
  471. C IF(DEP0+XINCR.GT.DPLUS0+OPLUS0)THEN
  472. C 2003 IF(DEP0+XINCR.GT.(UN-XPETIT)*DPLUS0)THEN
  473. IF(DEP0+XINCR.GT.(UN-XPETIT)*(DPLUS0+OPLUS0))THEN
  474. LSWICH=.TRUE.
  475. C 2003 XINCR =DPLUS0-DEP0
  476. XINCR =DPLUS0+OPLUS0-DEP0
  477. ENDIF
  478. C-->On incremente explicitement tout ce qui doit etre incremente
  479. C WARNING:l'increment de d.p. est calcule en 2 steps
  480. C WARNING:force effective pointee=FPLUS0*(UN-APLUS0)
  481. C 2003 FORF=FOR0+XINCR/(DPLUS0-DEP0)*(FPLUS0*(UN-APLUS0)-FOR0)
  482. FORF=FOR0+XINCR/(DPLUS0+OPLUS0-DEP0)*(FPLUS0*(UN-APLUS0)-FOR0)
  483. DEP0=DEP0+XINCR
  484. DENE=ABS(FCAMP0/2*(XINCR-(FORF-FOR0)/KCINE0))
  485. DDINL=DDINL+XINCR+FOR0/KCINE0
  486. FOR0=FORF
  487. C
  488. EPLUS0=EPLUS0+ ALPHP *DENE
  489. C 2003 EMOIN0=EMOIN0+(UN-ALPHP)*DENE
  490. EMOIN0=EMOIN0+ ALPHM *DENE
  491. * APLUS0=MIN(UN,APLUS0+BETA* ALPHP* DENE)
  492. * AMOIN0=MIN(UN,AMOIN0+BETA*(UN-ALPHP)*DENE)
  493. C 2003 APLUS0=TANH(BETA*EPLUS0)
  494. C 2003 AMOIN0=TANH(BETA*EMOIN0)
  495. CALL MSHEA3(BETA,EXPN, EPLUS0,EMOIN0, APLUS0,AMOIN0)
  496. C-->On calcule les caracteristiques de la surface
  497. C Cinematique en tenant compte de la zone de transition
  498. CALL MSHEA1(DEP0,ICAS,FOR0,DPLUS0,DMOIN0,
  499. > DTRANP,DTRANM,
  500. > APLUS0,CURFP,NCURFP,CURKP,NCURKP,CURLP,NCURLP,
  501. > AMOIN0,CURFM,NCURFM,CURKM,NCURKM,CURLM,NCURLM,
  502. > FCINI0,FCAMP0,KCINE0)
  503. C2003> FCINI0,FCAMP0,KCINE0,ALPHP)
  504. C+2003
  505. CC CALL MSHEA2(DEP0,DTRANP+OPLUS0,DTRANM+OMOIN0,
  506. CALL MSHEA2(DEP0,DTRANP+OPLUS0+OMOIN0,DTRANM+OPLUS0+OMOIN0,
  507. > FOR0,FCAMP0,KCINE0,ICAS,
  508. > ALPHP,ALPHM)
  509. C+2003
  510. DDINL=DDINL-FOR0/KCINE0
  511. C-->Switch sur B0
  512. IF(LSWICH)THEN
  513. ICAS=1
  514. VARF(3)=EPLUS0
  515. VARF(4)=EMOIN0
  516. VARF(6)=DMOIN0
  517. VARF(10)=APLUS0
  518. VARF(11)=AMOIN0
  519. GOTO 2
  520. ENDIF
  521. C-->Fin normale loop sur les sous-increments
  522. ENDDO
  523. C
  524. FORF=FOR0
  525. C
  526. VARF(1)=ICAS
  527. VARF(3)=EPLUS0
  528. VARF(4)=EMOIN0
  529. VARF(7)=FCINI0
  530. VARF(8)=FCAMP0
  531. VARF(9)=KCINE0
  532. VARF(10)=APLUS0
  533. VARF(11)=AMOIN0
  534. C
  535. RETURN
  536. C
  537. END
  538.  
  539.  
  540.  
  541.  

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