Télécharger sjonl3.eso

Retour à la liste

Numérotation des lignes :

sjonl3
  1. C SJONL3 SOURCE CHAT 06/06/01 21:20:08 5450
  2. SUBROUTINE SJONL3(sigi, depst, vari, xmat, tras, jtras ,
  3. > trat, jtrat, trac, jtrac, sigf, varf,
  4. > depsp, ireturn)
  5. c --- SJONL3 ---
  6. c computing plastic flow in joint element
  7. c under 3D plane (strain) conditions
  8. c
  9. c beat release 24-7-00 ylp.
  10. c sigi IN initial stress
  11. c depst IN total strain increment
  12. c vari IN initial internal variable
  13. c xmat IN material properties
  14. c tras IN shear plastic softening behaviour
  15. c jtras IN <tras> number of points
  16. c trat IN tensile plastic softening behaviour
  17. c jtrat IN <trat> number of points
  18. c trac IN compression plastic behaviour
  19. c jtrac IN <trac> number of points
  20. c sigf OUT final stress
  21. c varf OUT final internal variable
  22. c depsp OUT plastic strain increment
  23. c ireturn OUT routine control integer
  24.  
  25. IMPLICIT REAL*8(A-H,K-Z)
  26.  
  27. IMPLICIT INTEGER(I-J)
  28. dimension sigi(*), depst(*), vari(*),
  29. > sigf(*), depsp(*), varf(*)
  30. dimension ddepst(3), ddepsp(3)
  31. dimension xmat(*)
  32. dimension tras(2,jtras), trat(2,jtrat), trac(2,jtrac)
  33. c 1 --> strain
  34. c 2 --> stress
  35. c_______________________________________________________________________
  36. c
  37. * save iloop
  38. * data iloop/0/
  39. c print *,'debug -- sjonl3d -- check everything'
  40. c
  41. c do i=1,3
  42. c print *,'depst ',depst(i)
  43. c print *,'vari ', vari(i)
  44. c print *,'sigi ', sigi(i)
  45. c print *,'sigf ', sigf(i)
  46. c print *,'varf ',varf(i)
  47. c print *,'depst ',depst(i)
  48. c enddo
  49. c do i=1,6
  50. c print *,'xmat ', xmat(i)
  51. c enddo
  52. c print *,'tras ', tras ,' jtras ', jtras
  53. c print *,'trat ', trat ,' jtrat ', jtrat
  54. c print *,'trac ', trac ,' jtrac ', jtrac
  55. c
  56. c print *,'ireturn ', ireturn
  57. * iloop = iloop + 1
  58. c look at xmat(6) and deduce elastic/ plastic...
  59. c ...coupled/uncoupled variables.
  60. c all variables Jxxx are integer used as logical with the
  61. c following rule (0 is true, 1 is false)
  62. c JCP : use plastic law in compression
  63. c JSP : use plastic law in shear
  64. c JTP : use plastic law in tension
  65. c JST : couple shear and tension
  66. c JSC : couple shear and compression
  67. call sjoncpl(xmat(6),JCP,JSP,JTP,JST,JSC)
  68. * print *,'DEBUG --> entree dans sjonl2'
  69. * print *,'DepsN DepsS ' ,depst(2), depst(1)
  70. * set some usefull variables
  71. ZERO = 0.0d0
  72. ONE = 1.0d0
  73. c set the initial location of the edge of the cone
  74. ETi = xmat(5)
  75. c index definition
  76. Ipress = 2
  77. Ishear = 1
  78. Itotal = 3
  79.  
  80. Ishea1 = 1
  81. Ishea2 = 2
  82. Inorm = 3
  83. Icompr = 4
  84.  
  85. Istrain = 1
  86. Istress = 2
  87. c number of increments
  88. INCR = 5
  89. c______________________________________________________________________
  90. c
  91. c elastic initialisation of the joint stiffness --> K
  92. c K[z]
  93. c with z = t tensile ; s shear ; c compression
  94. Kt0 = xmat(Ipress)
  95. Ks0 = xmat(Ishear)
  96. Kc0 = Kt0
  97. c print *,'Deps - tension ',depst(Ipress)
  98. c print *,'Deps - shear ',depst(Ishear)
  99. c split the displacement increment in few smaller increments
  100. ddepst(Inorm) = depst(Inorm)/INCR
  101. ddepst(Ishea1) = depst(Ishea1)/INCR
  102. ddepst(Ishea2) = depst(Ishea2)/INCR
  103. c STARTING MAIN LOOP
  104. c [0
  105. do iflag =1 , INCR
  106. c______________________________________________________________________
  107. c
  108. c initial residual strength
  109. c RTS --> tensile
  110. c RSS --> shear
  111. c
  112. c internal variables initialisation --> Alpha
  113. c A[xyz]
  114. c with x = i initial ; f final ; k intermediate computational value
  115. c y = p plastic ; t total
  116. c z = t tensile ; s shear
  117. Aipt = vari(Ipress)
  118. c \--> initial plastic tensile strain
  119. Aips = vari(Ishear)
  120. c \--> initial plastic shear strain
  121. Aipc = vari(Icompr)
  122. c \--> initial plastic compression strain
  123. c______________________________________________________________________
  124. c
  125. c initial residual strength
  126. c RTS --> tensile
  127. c RSS --> shear
  128. c RCS --> compression
  129. c compute the tensile stress RTS corresponding to the initial
  130. c plastic tensile strain with <trat> evolution
  131. call yofxcu(Aipt ,trat, jtrat, RTS, dRTS, ireturn)
  132. c \--> Aipt initial plastic tensile strain
  133. c RTS OUT
  134. c dRTS = dTkt/dAkps useless value
  135. c [1
  136. if (ireturn .ne. 0) return
  137. c \--> ERROR
  138. c 1]
  139. c compute the initial cohesion residual
  140. call yofxcu(Aips, tras, jtras, RSS, dRSS, ireturn)
  141. c [1
  142. if (ireturn .ne. 0) return
  143. c \--> ERROR
  144. c 1]
  145. c comput the initial compression residual
  146. call yofxcu(Aipc, trac, jtrac, RCS, dRCS, ireturn)
  147. c [1
  148. if (ireturn .ne. 0) return
  149. c \--> ERROR
  150. c 1]
  151. c______________________________________________________________________
  152. c
  153.  
  154. Aitt = vari(Itotal)
  155. c \--> initial total tensile strain
  156. c compute final total tensile strain...
  157. Aftt = Aitt + ddepst(Inorm)
  158. c ... and the initial tensile strength...
  159. TT0 = trat(Istress,1)
  160. c ... and the elastic compressive strength
  161. TC0 = trac(Istress,1)
  162. c print *,'TC0 ',TC0
  163. c check one or two things if entering sjonl2 for first time
  164. c [1
  165. * mis en commentaire par TC
  166. * if ( iloop .eq. 1) then
  167. c [2
  168. * if ( TT0 .eq. ZERO ) then
  169. * print *,'cut-off tensile strength is zero'
  170. * print *,'denied option'
  171. * return
  172. * endif
  173. c 2]
  174. c [2
  175. * if ( TT0 .gt. ETi ) then
  176. * print *,'WARNING cut-off ouside the friction cone...'
  177. * print *,'... tensile strength at the edge of the cone'
  178. * print *,'set to cut-off tensile strength'
  179. * return
  180. * endif
  181. c 2]
  182. * endif
  183. * fin de mise en commentaire par TC
  184. c 1]
  185. c initialization of the damage variables
  186. Dtensi = 1.0d0
  187. Dshear = 1.0d0
  188. Dcompr = 1.0d0
  189. c WARNING the damage variables are understand here as :
  190. c D = 1 means no damage
  191. c D = 0 means full damage...
  192. c which is exactly the opposite of the usual definition !!!
  193. c compute the predicted elastic stress state
  194. c using the secant stiffness
  195. c
  196. c first, look at the tension modulus...
  197. c maximal (initial) tensile stress initialisation
  198. c N.B.: it is thus assumed that the behaviour in tension is elastic
  199. c plastic, and that the peak value features the end of the elastic
  200. c regime.
  201. c check whether tension or compression region is reached
  202. c [1
  203. if ( Aftt .gt. ZERO ) then
  204. c \--> tensile region reached...
  205. c [2
  206. if ( JTP .ne. 0 ) then
  207. c \--> use elastic law
  208. Kt = Kt0
  209. else
  210. c \--> use damaged elastic predictor
  211. c [3
  212. if ( RTS .eq. ZERO ) then
  213. c \--> zero residual tensile strength
  214. Kt = ZERO
  215. c \--> tensile modulus set to zero
  216. else
  217. c \--> still got residual tensile strength
  218. Kt = RTS / ( Aipt + TT0 / Kt0 )
  219. endif
  220. c 3]
  221. c update tensile damage variable to enable coupling effect
  222. c [3
  223. if ( JST .eq. 0) Dtensi = RTS / TT0
  224. c 3]
  225. endif
  226. c 2]
  227. c print *,'Dtensi ' ,Dtensi
  228. JCT = 0
  229. c JCT is an integer value used as a logical to control the location
  230. c of the stress state with respect to tension or compression:
  231. c tension --> JCT = 0
  232. c compression --> JCT = 1
  233. else
  234. c \--> compression region reached...
  235. c [2
  236. if ( JCP .ne. 0 ) then
  237. c \--> use elastic law
  238. Kc = Kc0
  239. else
  240. c \--> use damaged elastic predictor
  241. c [3
  242. if ( RCS .eq. ZERO ) then
  243. Kc = ZERO
  244. else
  245. Kc = RCS / (Aipc + TC0 / Kc0 )
  246. endif
  247. c 3]
  248. c update compression damage variable
  249. c Dcompr = Kc / Kc0 --> too fast
  250. c ... try something else...
  251. endif
  252. c 2]
  253. JCT = 1
  254. endif
  255. c ]1
  256. c ... then, look at the shear modulus
  257. TS0 = tras(Istress,1)
  258. c [1
  259. if (JSP .ne. 0 ) then
  260. Ks = Ks0
  261. else
  262. c [2
  263. if ( RSS .eq. ZERO ) then
  264. c \--> zero residual shear strength
  265. c or opened bonds
  266. Ks = ZERO
  267. c \--> shear modulus set to zero
  268. c Such case shall only appear if the bonds are opened.
  269. c otherwise, residual cohesion is generally observed
  270. else
  271. c \--> still got residual shear strength
  272. Ks = RSS / (Aips + TS0 / Ks0)
  273. endif
  274. c update the shear damage variable
  275. if (JST .eq. 0) Dshear = RSS / TS0
  276. c(Ks / Ks0)
  277. c print *,'Dshear ' , Dshear
  278. c 2]
  279. endif
  280. c 1]
  281.  
  282. c [1
  283. if (JCT .eq. 0) then
  284. c check if previous stress state was in compression
  285. c [2
  286. if ( (Aitt * Aftt) .lt. ZERO ) then
  287. Tft = Kt * Aftt * Dtensi
  288. else
  289. Tft = sigi(Inorm) + ( ddepst(Inorm) * Kt * Dshear)
  290. endif
  291. c 2]
  292. Tfs1 = sigi(Ishea1) + ( ddepst(Ishea1) * Ks * Dtensi)
  293. Tfs2 = sigi(Ishea2) + ( ddepst(Ishea2) * Ks * Dtensi)
  294. Tfs0 = sqrt((Tfs1**2)+(Tfs2**2))
  295. else
  296. c check if previous stress state was in tension
  297. c [2
  298. if ( (Aitt * Aftt) .lt. ZERO ) then
  299. Tfc = Kc * Aftt * Dcompr
  300. else
  301. c Tfc = sigi(Ipress) + ( ddepst(Ipress) * Kc * Dshear)
  302. Tfc = sigi(Inorm) + ( ddepst(Inorm) * Kc )
  303. endif
  304. c 2]
  305. Tfs1 = sigi(Ishea1) + ( ddepst(Ishea1) * Ks * Dcompr)
  306. Tfs2 = sigi(Ishea2) + ( ddepst(Ishea2) * Ks * Dcompr)
  307. Tfs0 = sqrt((Tfs1**2)+(Tfs2**2))
  308. endif
  309. c 1]
  310.  
  311. c print *,'stress state elastic prediction'
  312. c print *,'tension compression shear'
  313. c print *,Tft, Tfc, Tfs1, Tfs2
  314. c______________________________________________________________________
  315. c
  316. c before checking if any kind of yielding occured, update the
  317. c final internal variables in the case of simple elastic regime.
  318. Afpt = Aipt
  319. Afps = Aips
  320. Afpc = Aipc
  321. c______________________________________________________________________
  322. c
  323. c now, check yielding
  324. c N.B.: note that each mode might be activated independently
  325. c without any restriction!!!!
  326. c first, look at the tension/compression mode...
  327. c [1
  328. if ( JCT .eq. 0 ) then
  329. c [2
  330. if ( JTP .eq. 0 ) then
  331. c [3
  332. if ( abs(Tft) .gt. RTS*Dshear ) then
  333. c \--> YIELDING IN TENSION
  334. c print *,'yielding in tension...'
  335. c Afpt = Aftt - TT0 / Kt0
  336. Afpt = Aipt + (Tft-RTS*Dshear)/Kt
  337. c print *,'Afpt ',Afpt
  338. c \--> updated plastic tensile strain
  339. * print *,'updated plastic tensile strain ', Afpt
  340. c compute the updated residual tensile strength...
  341. call yofxcu(Afpt, trat, jtrat, RTS ,dRTS, ireturn)
  342. c \--> Afpt updated plastic tensile strain
  343. c RTS updated residual tensile strength
  344. c dTkt dummy
  345. Tft = RTS * Dshear
  346. c \--> compute the final tensile stress
  347. c print *,'residual tensile strength ',Tft
  348. endif
  349. c 3]
  350. endif
  351. c 2]
  352. else
  353. c [2
  354. if ( JCP .eq. 0 ) then
  355. c if ( abs(Tfc) .gt. abs(RCS)*Dshear) then
  356. c [3
  357. if ( abs(Tfc) .gt. abs(RCS) ) then
  358. c \--> YIELDING IN COMPRESSION
  359. c print *,'yielding in compression'
  360. c Afpc = abs(Aftt) - TC0 / Kc0
  361. c Afpc = Aipc + ((Tfc*-1.0d0)-RCS*Dshear)/Kc
  362. Afpc = Aipc + ((Tfc*(-1.0d0))-RCS)/Kc
  363. c \--> updated plastic compression strain
  364. c print *,'updated plastic compression strain ', Afpc
  365. c compute the updated residual compression strength
  366. call yofxcu(Afpc, trac, jtrac, RCS , dRCS, ireturn)
  367. c then, go back to the compression domain (negative stress)
  368. Tfc = RCS * (-1.0d0)
  369. c print *,'residual compressive strength ',RCS
  370. endif
  371. c 3]
  372. endif
  373. c 2]
  374. endif
  375. c 1]
  376. c ... finally, look at the shear mode...
  377. c compute the tensile value at the edge using an homothetic transform.
  378. c the cut-off value plays the role of a drag-stress
  379. c ETf = RTS + ETi - TT0
  380. c ... definitely not relevant!!!!
  381. c ... try something else... like a simple homothetic rule related to
  382. c the location of the residual cut-off with regard to the initial
  383. c location of the edge of the friction cone
  384. ETf = ETi * RTS / TT0
  385. c the behaviour in shear is assumed to follow a
  386. c homothetic rule with respect to the behaviour in pure shear
  387. c the coefficient of this homothetic transform is given by...
  388. c [1
  389. if ( JSP .eq. 0 ) then
  390. c [2
  391. if ( ETf .ne. ZERO ) then
  392. c [3
  393. if ( JCT .eq. 0 ) then
  394. c \--> shear yielding with tension
  395. Chom = (ETf - Tft)/(ETf)*Dtensi
  396. c print *,'Coef. hom. ', Chom
  397. c WARNING the original version used
  398. c Chom = (RTS - Tft)/(RTS +Pnor)
  399. c xmat(5) or PNOR is a simple numerical value used is order to
  400. c avoid division by zero if the residual cohesion vanishes.
  401. c print *,'???', Dtensi, ONE
  402. c [4
  403. if ( Dtensi .eq. ONE) then
  404. c \--> no tensile damage
  405. c [5
  406. if ( abs(Tfs0) .gt. (RSS*Chom) ) then
  407. c \--> YIELDING IN SHEAR
  408. c print *,'YIELDING IN SHEAR'
  409. c compute the total plastic shear strain...
  410. Afps = Aips + ((abs(Tfs0)/Chom)-RSS)/Ks
  411. c print *,'plasti shear strain ', Afps
  412. c ...and the final cohesion residual
  413. call yofxcu(Afps,tras,jtras,RSS,dRSS,ireturn)
  414. c [6
  415. if (ireturn .ne. 0) return
  416. c \--> ERROR
  417. c 6]
  418. Tfs = sign(Chom*RSS , Tfs)
  419. if (Tfs0 .eq. ZERO) then
  420. Tfs1 = ZERO
  421. Tfs2 = ZERO
  422. else
  423. Tfs1= Tfs*Tfs1/Tfs0
  424. Tfs2= Tfs*Tfs2/Tfs0
  425. endif
  426. c print *,'updated shear stress ',Tfs
  427. c print *,'Tfs1 Tfs2 ',Tfs1,Tfs2
  428. endif
  429. c 5]
  430. else
  431. c \--> already got some damage due to tension
  432. Afpstmp = Aips+ (abs(Tfs0)/Chom-RSS)/(Ks*Dtensi)
  433. c print *, Afpstmp
  434. c [5
  435. if ( Afpstmp .gt. ZERO ) then
  436. call yofxcu(Afpstmp,tras,jtras,RSSTMP,dRSS,ireturn)
  437. c [6
  438. if (ireturn .ne. 0) return
  439. c \--> ERROR
  440. c 6]
  441. c print *,'RSSTMP', RSSTMP
  442. c [6
  443. if ( abs(Tfs0) .gt. (Chom*RSSTMP) ) then
  444. c \--> YIELDING IN SHEAR
  445. Afps = Afpstmp
  446. print *,'YIELDING IN SHEAR 2'
  447. RSS = RSSTMP
  448. Tfs = sign(Chom*RSS , Tfs)
  449. if (Tfs0 .eq. ZERO) then
  450. Tfs1 = ZERO
  451. Tfs2 = ZERO
  452. else
  453. Tfs1= Tfs*Tfs1/Tfs0
  454. Tfs2= Tfs*Tfs2/Tfs0
  455. endif
  456. print *,'updated shear stress ',Tfs
  457. print *,'Tfs1 Tfs2 ',Tfs1,Tfs2
  458. endif
  459. c 6]
  460. endif
  461. c 5]
  462. endif
  463. c 4]
  464. else
  465. c \--> check shear yielding in compression region
  466. Chom = (ETf - Tfc)/(ETf)*Dcompr
  467. c print *,'Coef. hom. ', Chom
  468. c print *,'???', Dcompr, ONE
  469. c [4
  470. if ( Dcompr .eq. ONE) then
  471. c \--> no crushing damage
  472. c [5
  473. if ( abs(Tfs0) .gt. (RSS*Chom) ) then
  474. c \--> YIELDING IN SHEAR
  475. c print *,'YIELDING IN SHEAR'
  476. c compute the total plastic shear strain...
  477. Afps = Aips + ((abs(Tfs0)/Chom)-RSS)/Ks
  478. c print *,'plasti shear strain ', Afps
  479. c ...and the final cohesion residual
  480. call yofxcu(Afps,tras,jtras,RSS,dRSS,ireturn)
  481. c [6
  482. if (ireturn .ne. 0) return
  483. c \--> ERROR
  484. c 6]
  485. Tfs = sign(Chom*RSS , Tfs)
  486. if (Tfs0 .eq. ZERO) then
  487. Tfs1 = ZERO
  488. Tfs2 = ZERO
  489. else
  490. Tfs1= Tfs*Tfs1/Tfs0
  491. Tfs2= Tfs*Tfs2/Tfs0
  492. endif
  493. c print *,'updated shear stress ',Tfs
  494. c print *,'Tfs1 Tfs2 ',Tfs1,Tfs2
  495. endif
  496. c 5]
  497. else
  498. c \--> already got some damage due to compression
  499. Afpstmp = Aips+ (abs(Tfs0)/Chom-RSS)/(Ks*Dcompr)
  500. c print *, Afpstmp
  501. c [5
  502. if ( Afpstmp .gt. ZERO ) then
  503. call yofxcu(Afpstmp,tras,jtras,RSSTMP,dRSS,ireturn)
  504. c [6
  505. if (ireturn .ne. 0) return
  506. c \--> ERROR
  507. c 6]
  508. c print *,'RSSTMP', RSSTMP
  509. c [6
  510. if ( abs(Tfs0) .gt. (Chom*RSSTMP) ) then
  511. c \--> YIELDING IN SHEAR
  512. Afps = Afpstmp
  513. c print *,'YIELDING IN SHEAR 2'
  514. RSS = RSSTMP
  515. Tfs = sign(Chom*RSS , Tfs)
  516. if (Tfs0 .eq. ZERO) then
  517. Tfs1 = ZERO
  518. Tfs2 = ZERO
  519. else
  520. Tfs1= Tfs*Tfs1/Tfs0
  521. Tfs2= Tfs*Tfs2/Tfs0
  522. endif
  523. c print *,'updated shear stress ',Tfs
  524. c print *,'Tfs1 Tfs2 ',Tfs1,Tfs2
  525. endif
  526. c 6]
  527. endif
  528. c 5]
  529. endif
  530. c 4]
  531. endif
  532. c 3]
  533. else
  534. c \--> the edge has reach the origin... no homothetic rule
  535. c can be computed !
  536. print *,'EDGE ZERO -- ERROR'
  537. endif
  538. c 2]
  539. endif
  540. c 1]
  541. c______________________________________________________________________
  542. c
  543. c final stage --> updating the output variables...
  544. c ... first the final stress...
  545. c [1
  546. if ( JCT .eq. 0 ) then
  547. sigf(Inorm) = Tft
  548. sigi(Inorm) = Tft
  549. else
  550. sigf(Inorm) = Tfc
  551. sigi(Inorm) = Tfc
  552. endif
  553. c 1]
  554. sigf(Ishea1) = Tfs1
  555. sigi(Ishea1) = Tfs1
  556. sigf(Ishea2) = Tfs2
  557. sigi(Ishea2) = Tfs2
  558. c ... then the final internal variables...
  559. varf(Ipress) = Afpt
  560. varf(Ishear) = Afps
  561. varf(Icompr) = Afpc
  562. varf(Itotal) = Aftt
  563. c ... and finally, the plastic strain increments.
  564. c [1
  565. if ( JCT .eq. 0 ) then
  566. ddepsp(Inorm) = varf(Ipress) - vari(Ipress)
  567. else
  568. ddepsp(Inorm) = varf(Icompr) - vari(Icompr)
  569. endif
  570. c 1]
  571. ddepsp(Ishea1) = varf(Ishear) - vari(Ishear)
  572. ddepsp(Ishea2) = varf(Ishear) - vari(Ishear)
  573. vari(Ipress) = Afpt
  574. vari(Ishear) = Afps
  575. vari(Icompr) = Afpc
  576. vari(Itotal) = Aftt
  577. depsp(Inorm) = depsp(Inorm) + ddepsp(Inorm)
  578. depsp(Ishea1) = depsp(Ishea1) + ddepsp(Ishea1)
  579. depsp(Ishea2) = depsp(Ishea2) + ddepsp(Ishea2)
  580. c END OF MAIN LOOP
  581. c 0]
  582. enddo
  583. c print *,'output variables :'
  584. c print *,' --> stress state - normal shear'
  585. c print *, sigf(Inorm), sigf(Ishea1), sigf(Ishea2)
  586. c print *,'internal variables - normal - comp - shear - total'
  587. c print *, varf(Ipress), varf(Icompr),varf(Ishear), varf(Itotal)
  588. c print *,'plastic strain - normal - shear'
  589. c print *, depsp(Inorm), depsp(Ishea1),depsp(Ishea2)
  590. c And, that's all folks !!!!
  591. c print *,'sjonl2 ended normaly...'
  592. end
  593.  
  594.  
  595.  
  596.  

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